Dual-slope Method for Enhanced Depth Sensitivity in Diffuse Optical Spectroscopy

Abstract

An apparatus for earning out near-infrared spectroscopy using intensity-modulated near-infrared radiation or pulsed near-infrared radiation includes sources and detectors. For each source, there exists first and second distances. The first distance is a distance between the source and a first detector. The second distance is a distance between the source and the second detector. For each source, the difference between these two distances is the same. Additionally, wherein, for each source, the detector at a shorter distance is the same detector that is at a longer distance for the other source. A processor derives, from signals received by the detectors, a parameter indicative of two matched slopes. Tins parameter is either phase of the intensity-modulated near-infrared radiation or mean time-of-flight data for the pulsed near-infrared radiation. The processor then provides output data based on an average of the matched slopes. This promotes reduced sensitivity to superficial layers and enhanced sensitivity to deeper portions of a medium that is under investigation.

Claims

1. An apparatus comprising a processor and a spectrometer, said spectrometer comprising a first source, a second source, a first detector, and a second detector, wherein said processor is configured to receive data indicative of signals detected by said detectors in response to illumination by said sources, wherein said sources emit near-infrared radiation, wherein, for each source, there exists a first distance and a second distance, wherein said first distance is a distance between said source and said first detector, wherein said second distance is a distance between said source and said second detector, wherein there exists a difference between said first and second distances, wherein, for each source, said difference is the same, wherein, for each source, the detector at a shorter distance is the same detector that is at a longer distance for the other source, wherein said processor is configured to derive, from said signals, a parameter indicative of two matched slopes, wherein said near-infrared radiation is selected from the group consisting of intensity-modulated near-infrared radiation and pulsed near-infrared radiation, wherein said parameter is selected from the group consisting of phase of said near-infrared radiation and mean time-of-flight data for said near-infrared radiation, and wherein said processor is further configured to provide output data based on an average of said matched slopes, thereby promoting reduced sensitivity to superficial layers and increased sensitivity to deeper portions of a medium that is under investigation.

2. The apparatus of claim 1, wherein said sources and detectors are arranged along a linear array.

3. The apparatus of claim 1, wherein said sources and detectors are arranged at vertices of a parallelogram.

4. The apparatus of claim 1, wherein said sources and detectors are arranged at vertices of a trapezoid.

5. The apparatus of claim 1, wherein said sources and detectors are arranged at vertices of a three-segment line.

6. The apparatus of claim 1, wherein said sources and detectors are arranged at vertices of a rectangle.

7. The apparatus of claim 1, further comprising additional detectors, wherein said first and second detectors and said additional detectors are arranged in an array.

8. The apparatus of claim 1, further comprising a cuff for applying a periodic pressure.

9. The apparatus of claim 1, wherein said processor is configured to carry out coherent hemodynamics spectroscopy.

10. The apparatus of claim 1, wherein said near-infrared radiation is intensity-modulated near-infrared radiation.

11. The apparatus of claim 1, wherein said near-infrared radiation is pulsed near-infrared radiation.

12. The apparatus of claim 1, wherein said near-infrared radiation is intensity-modulated near-infrared radiation and wherein said parameter is selected to be a phase of said intensity-modulated near-infrared radiation.

13. The apparatus of claim 1, wherein said parameter is selected to be mean time-of-flight for said pulsed near-infrared radiation.

14. The apparatus of claim 1, wherein the number of detectors is odd.

15. The apparatus of claim 1, wherein the number of detectors is even.

16. A method comprising illuminating a region of skin with near-infrared radiation emitted by first and second source, detecting radiation with first and second detectors, wherein said sources emit intensity-modulated near-infrared radiation or pulsed near-infrared radiation, wherein, for each source, there exists a first distance, which is a distance between said source and said first detector, and a second distance, which is a distance between said source and said second detector, wherein there exists a difference between said first and second distances, wherein, for each source, said difference is the same, wherein, for each source, the detector at a shorter distance is the same detector that is at a longer distance for the other source, deriving, from data indicative of said radiation, a parameter indicative of two matched slopes, said parameter being one of phase of said intensity-modulated near-infrared radiation and mean time-of-flight data for said pulsed near-infrared radiation, and providing output data based on an average of said matched slopes, thereby promoting reduced sensitivity to superficial layers and increased sensitivity to a deeper portion of a medium that is under investigation.

17. The method of claim 16, further comprising selecting said region to be a scalp and selecting said deeper portion is a cerebral cortex.

18. The method of claim 16, further comprising selecting said parameter to be said phase of said intensity-modulated near-infrared radiation.

19. The method of claim 16, further comprising selecting said parameter to be said mean time-of-flight data for said pulsed infrared radiation.

20. An apparatus comprising a processor and a spectrometer, said spectrometer comprising one source and one detector, wherein said processor is configured to receive data indicative of signals detected by said detector in response to illumination by said source, wherein said source emits near-infrared radiation, wherein said near-infrared radiation is selected from the group consisting of intensity-modulated near-infrared radiation and pulsed near-infrared radiation, wherein said detector detects a parameter selected from the group consisting of phase of said intensity-modulated near-infrared radiation and mean time-of-flight of said pulsed near-infrared radiation, and wherein said processor is configured to derive, from said parameter, data indicative of a change in absorption properties of a medium under investigation, thereby promoting reduced sensitivity to superficial layers and increased sensitivity to deeper portions of said medium.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0030] FIG. 1A shows a system for carrying out a dual-slope method as described herein.

[0031] FIG. 1B shows one embodiment of the spectrometer shown in FIG. 1A;

[0032] FIG. 1C shows the spectrometer of FIG. 1C with exemplary spacings shown;

[0033] FIG. 1D shows an alternative configuration of sources and detectors for the spectrometer shown in FIG. 1A;

[0034] FIG. 2 shows sensitivity for layered changes in the absorption coefficient;

[0035] FIGS. 3 and 4 compares ratios of phasors;

[0036] FIG. 5 is a sensitivity map;

[0037] FIG. 6 is a section through the sensitivity map of FIG. 5;

[0038] FIG. 7 is another sensitivity map;

[0039] FIG. 8 is a section through the sensitivity map of FIG. 7;

[0040] FIG. 9 shows sensitivity maps for dual-slope data for intensity and phase;

[0041] FIG. 10 compares cross sections of sensitivity maps of FIGS. 7 and 9;

[0042] FIG. 11 compares data in FIGS. 7 and 8;

[0043] FIG. 12 shows sensitivity maps for dual-slope data;

[0044] FIG. 13 shows a diffusive medium;

[0045] FIGS. 14-15 show estimated phasors.

[0046] FIGS. 16 and 17 shows maps of signal-to-noise ratios; and

[0047] FIG. 18 shows a table of signal-to-noise ratios.

DETAILED DESCRIPTION

[0048] FIG. 1A shows a subject 10 under examination by a spectrometer 12 that is being operated by a processor 14. Preferably, the spectrometer 12 is a frequency-domain near-infrared spectrometer 12.

[0049] The subject 10 is shown as interacting with additional devices, such as a blood-pressure monitor 16 and an inflatable cuff 18. Although these are optional components, data provided by them is potentially useful when inspected in light of data provided by the spectrometer 12. For example, using the inflatable cuff 18, it is possible to apply a periodic stimulation that will permit the use of coherent hemodynamics spectroscopy. These additional devices are likewise in communication with the processor 14.

[0050] Referring now to FIG. 1B, the spectrometer 12 includes a symmetrical source detector array 18 having first and second sources 20, 22 and first and second detectors 24, 26 arranged along a line.

[0051] Referring now to FIG. 1C, the distance between the first source 20 and the first detector 24 shall be referred to herein as a “short distance.” This short distance is equal to the distance between the second source 22 and the second detector 26. The distance between the first source 20 and the second detector 26, referred to herein as the “long distance,” is equal to that between the second source 22 and the first detector 24. In the illustrated embodiment, the short distance is equal to twenty millimeters and the long distance is equal to thirty-five millimeters.

[0052] In principle, what needs to be equal is the difference between the long and short distances associated with each source. Thus, it is also possible to implement the array 18 as a two-dimensional array in which sources 20, 22 and detectors 24, 26 are, instead, on vertices of a rectangle.

[0053] Another embodiment, shown in FIG. 1D, features first and second sources 20, 22 and multiple detectors 28, 30, 32, 34, 36, 38, 40 with the same separations between sources and detectors.

[0054] In either case, the array 18 collects data from which it is possible to identify two matched slopes, one of which is associated with the first source and the other of which is associated with the second source. In the case of matched slopes, the roles of each source 20, 22 and each detector 24, 26 are switched such that any given source 20, 22 or detector 24, 26 that participates in generating data at the short distance simultaneously participates in generating data at the long distance. The array 18 thus provides a basis for implementing a method that measures two slopes and averages them. This method is referred to herein as the “dual-slope method.”

[0055] This dual-slope method features: (1) insensitivity to instrumental effects related to temporal variations in source emission and detector sensitivity properties; (2) insensitivity to changes in optical coupling between optical probe and tissue; (3) reduced sensitivity to localized as well as uniform superficial tissue inhomogeneities, resulting in a greater relative sensitivity to deeper tissue; and (4) localized sensitivity to a deep tissue volume, which does not feature the typical banana shape of near-infrared spectroscopy regions-of-sensitivity.

[0056] Disclosed herein is a detailed theoretical study of the sensitivity of single-distance vs. single-slope and dual-slope methods to focal absorption perturbations consistent with promoting achievement of points (3) and (4) above. The procedure is described as being carried out using any of three data types: DC, AC, and phase. The results suggest that the dual-slope method provides an effective way to obtain a signal from deep tissue with minimal sensitivity to superficial tissue and instrumental artifacts. In particular, the dual-slope method provides a way to obtain reliable data concerning the cortex of the brain without that data being significantly confounded by phenomena occurring in the scalp.

[0057] In particular, the results shown herein suggest that phase dual-slopes feature a more specific sensitivity to deeper tissue than intensity dual-slopes. However, a drawback of phase measurements is their lower signal-to-noise ratio compared to intensity measurements. For cases of lower optical contrast, or for instruments that do not collect frequency-domain data, dual DC intensity slopes may be used as they are still preferable over single-slopes, but they would feature maximal sensitivity to shallower tissue than dual phase slopes. Other theoretical computations were run to establish if the dual-slope method can correctly retrieve relative changes of oxy- and deoxyhemoglobin concentrations occurring in the brain in typical brain studies. Two scenarios are discussed: (a) a typical protocol of coherent hemodynamic spectroscopy and (b) a typical brain activation case.

[0058] In the first case, systemic hemodynamic changes are induced both in the extracerebral layer and in the brain rather uniformly so that simulations rely on the assumption of layered changes in the absorption coefficients. The second case relies on an assumption of more localized absorption changes deeper in the tissue and some absorption changes in the outer layer. The results described herein suggest that the dual-slope method, in particular with phase data, can effectively retrieve the relative changes of oxy- and deoxyhemoglobin concentrations in deeper tissues.

2. Theory

[0059] Section 2.1 introduces formulas useful for understanding the sensitivity of DC (direct current intensity), AC (alternating current amplitude), and phase (φ), measured at a single source-detector separation, to focal absorption changes within an arbitrary diffusing and absorbing medium. Section 2.2. derives the sensitivity of the single-slopes of these data types to focal absorption changes where the independent variable is the source-detector separation.

[0060] The derived formulas are agnostic to the geometry and optical heterogeneity of the medium. In a geometry of semi-infinite homogeneous diffusive medium, the quantities ln(r.sup.2DC(r)), ln(r.sup.2AC(r)), and φ(r), where r is the source-detector separation, are well approximated by straight lines. An assumption is made that these properties are also valid for many types of real tissues if the measurements are carried out in a reflectance geometry. Such a geometry is relied upon herein unless otherwise stated. In general the slope of a straight line and the way it changes in dynamic and heterogenous conditions depends on both the source-detector separations and on the source-detector arrangement on the tissue surface.

[0061] Section 2.3 addresses the sensitivity of a data type measured at a single-distance or its slope to focal hemoglobin changes. Section 2.4 introduces the method of the dual-slope with two different types of source-detector arrangements. Section 2.5 provides the equations used for the numerical results of this work. Section 2.6 discusses noise-to-signal ratio for the different methods.

[0062] In the following sections, a variable in bold means that the variable represents a position vector in the space, while an arrow on top of a variable indicates that the variable is not a vector in space but an array of a specified dimension.

2.1 Spatial Sensitivity of Raw Data Types (DC, AC, and Phase) at a Single Source-Detector Separation

[0063] Within first-order perturbation theory of the diffusion equation, one can write:

[00001]$\begin{array}{cc}\frac{1}{DC\left(r_s,r_d\right)}\frac{\partial DC\left(r_s,r_d\right)}{\partial {\mu }_{\mathrm{ak}}}=-.Math.l.Math.\left(r_s,r_k,r_d\right)=-.Math.l_{DC}.Math.\left(r_s,r_k,r_d\right)& \left(1\right)\end{array}$

where r.sub.s is the position vector of a point source, such as that associated with the first or second source 20, 22 and r.sub.d is the position vector associated with a detector 24, 26, DC (r.sub.s, r.sub.d) is the detected direct current intensity at a detector position r.sub.d, μ.sub.ak is the absorption coefficient of a region (labeled k) within the diffusive medium which is identified by the position vector r.sub.k, and l(r.sub.s, r.sub.k, r.sub.d) is the average partial pathlength traveled inside region k by photons emitted at r.sub.s by a particular source 20, 22 detected at r.sub.d by a particular detector 24, 26. It is possible to reframe the perturbation theory of the diffusion equation by using the pathlength moments in the three domains, namely the continuous-wave domain, the frequency domain, and the time domain, of near-infrared spectroscopy and to show that a formula similar to equation (1) is valid also for the frequency-domain reflectance:

[00002]$\begin{array}{cc}\frac{1}{\tilde{R}\left(r_s,r_d,\omega \right)}\frac{\partial \tilde{R}\left(r_s,r_d,\omega \right)}{\partial {\mu }_{\mathrm{ak}}}=-.Math.\tilde{l}.Math.\left(r_s,r_k,r_d,\omega \right)& \left(2\right)\end{array}$

where {tilde over (R)}(r.sub.s, r.sub.d, ω)=AC(r.sub.s, r.sub.d, ω)exp(r.sub.s, r.sub.d, ω)) is the complex notation for the reflectance measured by a detector 24, 26 at r.sub.d, and l is the associated complex pathlength, which reduces to the usual pathlength of equation (1) when the angular modulation frequency (ω) is zero. In fact, equation (2) reduces to equation (1) for ω=0. Equation (2) implies two separate equations for AC and φ:

[00003]$\begin{array}{cc}\frac{1}{AC\left(r_s,r_d,\omega \right)}\frac{\partial AC\left(r_s,r_d,\omega \right)}{\partial {\mu }_{\mathrm{ak}}}=-\mathrm{Re}\left[.Math.\tilde{l}.Math.\left(r_s,r_k,r_d,\omega \right)\right]=-.Math.l_{AC}.Math.\left(r_s,r_k,r_d,\omega \right)& \left(3\right)\\ \frac{\partial \phi \left(r_s,r_d,\omega \right)}{\partial {\mu }_{\mathrm{ak}}}=-\mathrm{Im}\left[.Math.\tilde{l}.Math.\left(r_s,r_k,r_d,\omega \right)\right]=-.Math.l_{\phi }.Math.\left(r_s,r_k,r_d,\omega \right)& \left(4\right)\end{array}$

An approximate formula for {tilde over (l)}(r.sub.s, r.sub.k, r.sub.d, ω) for the foregoing is:

[00004]$\begin{array}{cc}.Math.\tilde{l}.Math.\left(r_s,r_k,r_d,\omega \right)\cong \frac{\tilde{\Phi }\left(r_s,r_k,\omega \right)\tilde{R}\left(r_k,r_d,\omega \right)}{\tilde{R}\left(r_s,r_d,\omega \right)}{V}_k& \left(5\right)\end{array}$

where {tilde over (Φ)}(r.sub.s, r.sub.k, ω), {tilde over (R)}(r.sub.k, r.sub.d, ω), and {tilde over (R)}(r.sub.s, r.sub.d, ω) are the frequency-domain Green's function of fluence calculated at r.sub.k when photons are emitted from r.sub.s, the frequency-domain reflectance (output intensity) measured at the detector's point r.sub.d when photons are emitted from r.sub.k, and the frequency-domain reflectance measured at the detector's point r.sub.d when photons are emitted from r.sub.s, respectively. V.sub.k is the volume of the k-region (“centered” at r.sub.k) at which a change in absorption occurs. Equation (5) is valid if the maximum linear size of the k-region is much smaller than |r.sub.s−r.sub.k| and |r.sub.d−r.sub.k|. If this approximation is not valid, the correct calculation of the partial pathlength can instead by carried out by integrating over the relevant volume.

[0064] The complex partial pathlength can be interpreted as the fraction of the photon density wave detected at the detector 24, 26 at r.sub.d that reached the region at r.sub.k. Because photon density waves are exponentially damped as a function of the modulation frequency, increasing ω causes a smaller wave fraction to reach the k-region before being detected. On the contrary, for the limiting case ω=0 (DC case) the partial pathlength depends only on geometrical conditions and optical properties. It is useful to rewrite equations (1), (3), and (4) in a concise way with a single formula, in which the dependence on ω is implicit:

[00005]$\begin{array}{cc}\frac{\stackrel{\_}{\partial Y}\left(r_s,r_d\right)}{\partial {\mu }_{\mathrm{ak}}}=-.Math.l_Y.Math.\left(r_s,r_k,r_d\right)& \left(6\right)\end{array}$

where “Y” is a data type (DC, AC, or φ), the bar on the top of the derivative means that the derivative is normalized, in the DC and AC cases, by dividing by DC or AC, respectively, as in equations (1) and (3), and l.sub.Y(r.sub.s, r.sub.k, r.sub.d) is the pathlength corresponding to the Y data type, as in equations (1), (3), and (4), when the focal change occurs at r.sub.k, and the photons are emitted at r.sub.s and detected at r.sub.d.

[0065] For M distinct focal changes in the absorption coefficient, within first order perturbation theory (which assumes independent, i.e. non-interacting focal absorption changes) it is possible to express the total change of the normalized data type as the sum of contributions due to each focal change:

[00006]$\begin{array}{cc}\stackrel{\_}{\Delta Y}\left(r_s,r_d\right)=-\underset{k=1}{\overset{M}{.Math.}}.Math.l_Y.Math.\left(r_s,r_k,r_d\right)\Delta {\mu }_{\mathrm{ak}}& \left(7\right)\end{array}$

[0066] Many sophisticated algorithms target the reconstruction of the true focal absorption changes Δμ.sub.ak by measuring ΔY(r.sub.s, r.sub.d) for multiple source-detector pairs on the medium boundary and using computationally complex inversion procedures. The disclosure considers the “equivalent” homogeneous absorption change that yields the same ΔY(r.sub.s, r.sub.d) that is associated with a collection of M localized absorbers. In fact, based on equation (6), if one assumes that ΔY(r.sub.s, r.sub.d) arises from a homogeneous absorption change in the medium:

[00007]$\begin{array}{cc}\stackrel{\_}{\Delta Y}\left(r_s,r_d\right)=-.Math.L_Y.Math.\left(r_s,r_d\right)\Delta {\mu }_{aY}\left(r_s,r_d\right)& \left(8\right)\end{array}$

where L.sub.Y(r.sub.s, r.sub.d) is the average total pathlength travelled (in the whole medium) by photons emitted at r.sub.s and detected at r.sub.d for the data type Y, and Δμ.sub.aY(r.sub.s, r.sub.d) is the equivalent homogeneous absorption change estimated by the data type Y (when Y=DC, equation (8) is the modified Beer-Lambert law). The equivalent homogeneous absorption change is related to the set of M true absorption changes by the relationship:

[00008]$\begin{array}{cc}\Delta {\mu }_{\mathrm{aY}}\left(r_s,r_d\right)=\underset{k=1}{\overset{M}{.Math.}}\frac{.Math.l_Y.Math.\left(r_s,r_k,r_d\right)}{.Math.L_Y.Math.\left(r_s,r_d\right)}{\mathrm{\Delta \mu }}_{\mathrm{ak}}& \left(9\right)\end{array}$

[0067] For a single focal change in absorption at r.sub.k, the sensitivity, S, of the data type Y detected at r.sub.d is as follows:

[00009]$\begin{array}{cc}{S}_Y\left(r_s,r_k,r_d\right)=\frac{.Math.l_Y.Math.\left(r_s,r_k,r_d\right)}{.Math.L_Y.Math.\left(r_s,r_d\right)}=\frac{\Delta {\mu }_{\mathrm{aY}}\left(r_s,r_d\right)}{\Delta {\mu }_{\mathrm{ak}}}& \left(10\right)\end{array}$

[0068] S.sub.Y(r.sub.s, r.sub.k, r.sub.d) is always positive for DC data because l.sub.DC(r.sub.s, r.sub.k, r.sub.d) and L.sub.DC(r.sub.s, r.sub.d) represent the actual physical mean partial and total pathlengths, respectively. This is true also for AC data, at least for typical values of co used in near-infrared spectroscopy, even though l.sub.AC(r.sub.s, r.sub.k, r.sub.d) and L.sub.AC(r.sub.s, r.sub.d) do not represent mean physical pathlengths. But this is not generally true for phase data. Therefore, given a single focal change (μμ.sub.ak), the equivalent homogeneous absorption change estimated with DC and AC data (Δμ.sub.aDC, Δμ.sub.aAC) will always have the same sign of the true change, while this is not always the case for phase data.

[0069] The behavior of phase data with respect to a focal absorption change can be understood if one considers that, for typical modulation frequencies used in near-infrared spectroscopy such that ω<<vμ.sub.a (where v is the speed of light in the medium), φ(r.sub.s,r.sub.d)≅(ω/v)L.sub.DC(r.sub.s,r.sub.d). If one considers a focal absorption increase occurring deeper in tissue, photons travelling along longer paths will have a higher chance of being absorbed. The overall effect will be a shortening of L.sub.DC(r.sub.s, r.sub.d) and therefore a decrease in phase. If an absorption increase occurs closer to the medium boundary, the effect on phase is more complex. In such cases, it is possible that photons with shorter pathlengths will be preferentially absorbed. This would increase L.sub.DC(r.sub.s, r.sub.d) and therefore increase phase.

[0070] A more rigorous derivation is based on the observation that the change in phase due to a focal change in absorption depends on both the variance of the pathlengths inside the focal region and the covariance of the pathlengths travelled inside and outside the region. While the variance is always positive, the covariance can be either positive or negative and may prevail over the other term. However, a homogenous increase in absorption will always cause a decrease in phase. These properties can be summarized by saying that, in the case of an absorption increase: l.sub.φ(r.sub.s, r.sub.k, r.sub.d) may be >0, <0, or =0, whereas L.sub.φ(r.sub.s, r.sub.d) is always >0. The fact that L.sub.φ is always >0 is linked to the definition of L.sub.DC, where the pathlength of each detected photon (L) is weighted by the factor e.sup.−μ.sup.a.sup.L, which penalizes photons travelling longer paths with respect to photons travelling shorter paths. Therefore, as one uniformly increases the absorption coefficient of the whole medium, L.sub.DC decreases, and so does the phase. This statement can be proven mathematically by using the properties of the radiative-transfer equation.

[0071] All the equations herein are general and do not depend on the geometry of the medium (except media with concave boundaries) and heterogeneity of the optical properties. The only assumption is the validity of the first order perturbation theory of the diffusion equation.

[0072] In a practical situation, the equivalent absorption change is found by inverting equation (8). This requires the knowledge of the total pathlengths. For particular geometries, like the semi-infinite homogeneous medium geometry, the total pathlengths can be calculated once one knows the absolute optical properties of the medium (see § 2.5).

2.2 Spatial Sensitivity of Single-Slopes of DC, AC, and Phase

[0073] This section simplifies the notation by using scalar source-detector distances (r.sub.i) instead of position vectors for one source (r.sub.s) and multiple detectors (r.sub.di where i=1, . . . , N), which are on the tissue surface. In other words, r.sub.i=|r.sub.s−r.sub.di| will be used instead of the whole set of position vectors for the source and the detectors. However, in a heterogeneous medium, as assumed herein, the partial and total pathlengths when photons are collected at r.sub.di depends on both r.sub.s and r.sub.di (and on r.sub.k for the partial pathlengths) and not only on the distance r.sub.i. The 3-dimensional (3D) position vector notation is retained to identify the focal regions where absorption changes occur.

[0074] The following discussion provides formulas for the sensitivity of the slopes of ln(r.sup.2DC(r)), ln(r.sup.2AC(r)), and φ(r) to focal absorption changes.

[0075] The general solution for the best straight line (in a least square sense) through N points having coordinates (x.sub.i, y.sub.i) where i=1, . . . , N yields a slope equal to cov({right arrow over (x)}, {right arrow over (y)})/var({right arrow over (x)}), where “cov” and “var” are the covariance and variance, respectively. In the embodiments described herein, x.sub.i=r.sub.i is a set of N source-detector separations (with respect to a single source) and y.sub.i=y(r.sub.i) is ln [r.sub.i.sup.2DC(r.sub.i)], ln [r.sub.i.sup.2AC(r.sub.i)], or φ(r.sub.i), which can be considered as new data types (“y”) derived from the data types “Y”. More specifically, this section considers slopes that are obtained using a single source and multiple detectors (or, equivalently, a single detector and multiple sources) and refers to these slopes as “single-slopes.” Therefore, by indicating the single-slope of data type “y” as SSl.sub.Y, one can write:

[00010]$\begin{array}{cc}SS{l}_Y=\frac{\mathrm{cov}\left(\stackrel{.fwdarw.}{r},\stackrel{.fwdarw.}{y}\right)}{\mathrm{var}\left(\stackrel{.fwdarw.}{r}\right)}& \left(11\right)\end{array}$

Differentiating equation (11) with respect to a focal change in absorption at (r.sub.k) yields:

[00011]$\begin{array}{cc}\frac{\partial {\mathrm{SSl}}_Y\left(r_k,\stackrel{.fwdarw.}{r}\right)}{\partial {\mu }_{\mathrm{ak}}}=-\frac{\left(\stackrel{.fwdarw.}{r}-\stackrel{.fwdarw.}{.Math.\stackrel{.fwdarw.}{r}.Math.}\right)\stackrel{.fwdarw.}{.Math.l_Y.Math.}\left(r_k,\stackrel{.fwdarw.}{r}\right)}{N\mathrm{var}\left(\stackrel{.fwdarw.}{r}\right)}& \left(12\right)\end{array}$

where {right arrow over (r)} is the average source-detector separation ({right arrow over (r)}=Σ.sub.i=1.sup.N r.sub.i/N), ({right arrow over (r)}−) is a 1×N array with elements r.sub.i−{right arrow over (r)}, (r.sub.k, {right arrow over (r)}) is an N×1 array with elements l.sub.Y(r.sub.k, r.sub.i), and var({right arrow over (r)})=(|{right arrow over (r)}−|.sup.2).

[0076] The expression l.sub.Y(r.sub.k, r.sub.i) is a simplified notation for l.sub.Y(r.sub.s, r.sub.k, r.sub.di) where i=1, . . . , N. The derivation of equation (12) uses equations (2)-(4) and the property that

[00012]$\frac{\partial \mathrm{yi}}{\partial {\mu }_{\mathrm{ak}}}=\left(\frac{\stackrel{\_}{\partial Y}\left(r_i\right)}{\partial {\mu }_{\mathrm{ak}}}\right)=-.Math.l_Y.Math.\left(r_k,r_i\right).$

Within first order perturbation theory, it is possible to derive the change in the slope of a data type due to M focal absorbers at r.sub.k where k=1, . . . , M:

[00013]$\begin{array}{cc}\Delta {\mathrm{SSl}}_Y\left(\stackrel{.fwdarw.}{r}\right)=-\frac{\left(\stackrel{.fwdarw.}{r}-\stackrel{.fwdarw.}{.Math.\stackrel{.fwdarw.}{r}.Math.}\right)\left(\right)\stackrel{.fwdarw.}{{\mathrm{\Delta \mu }}_a}}{N\mathrm{var}\left(\stackrel{.fwdarw.}{r}\right)}& \left(13\right)\end{array}$

where is an N×M matrix with elements l.sub.Y(r.sub.k, r.sub.i), and {right arrow over (Δμ.sub.a)} is an M×1 absorption change array (with elements Δμ.sub.ak).

[0077] The expression (r.sub.k, {right arrow over (r)}) of equation (12) is the k.sup.th column of . As was the case for equation (7) for the raw data types, equation (13) shows how the changes in the slope of a data type are related to the true localized changes in the absorption coefficients in the medium. Also, for the slopes, it is possible to estimate an equivalent homogeneous absorption change (Δμ.sub.aY) associated with ΔSSl.sub.Y which satisfies the equation:

[00014]$\begin{array}{cc}\Delta {\mathrm{SSl}}_Y\left(r_s,\stackrel{.fwdarw.}{r}\right)=-\frac{\left(\stackrel{.fwdarw.}{r}-\stackrel{.fwdarw.}{.Math.\stackrel{.fwdarw.}{r}.Math.}\right)\stackrel{.fwdarw.}{.Math.L_Y.Math.}\left(\stackrel{.fwdarw.}{r}\right)\Delta {\mu }_{aY}\left(\stackrel{.fwdarw.}{r}\right)}{N\mathrm{var}\left(\stackrel{.fwdarw.}{r}\right)}& \left(14\right)\end{array}$

where ({right arrow over (r)}) is an N×1 array with elements L.sub.Y(r.sub.i) where i=1, . . . , N. In this case, the homogeneous equivalent absorption change Δμ.sub.aY({right arrow over (r)}) depends not only on the data type (Y) but also on the set of source-detector distances ({right arrow over (r)}). The relationship between the true focal absorption changes and Δμ.sub.aY is easily derived as:

[00015]$\begin{array}{cc}\Delta {\mu }_{\mathrm{aY}}\left(\stackrel{.fwdarw.}{r}\right)=\frac{\left(\stackrel{.fwdarw.}{r}-\stackrel{.fwdarw.}{.Math.\stackrel{.fwdarw.}{r}.Math.}\right).Math..Math.\stackrel{.fwdarw.}{{\mathrm{\Delta \mu }}_a}}{\left(\stackrel{.fwdarw.}{r}-\stackrel{.fwdarw.}{.Math.\stackrel{.fwdarw.}{r}.Math.}\right)\stackrel{.fwdarw.}{.Math..Math.}\left(\stackrel{.fwdarw.}{r}\right)}& \left(15\right)\end{array}$

Based on equation (15), it is possible to define the sensitivity of the single-slope of a data type (where the set of source-detector distances {right arrow over (r)} is used) with respect to a focal absorption change at r.sub.k as:

[00016]$\begin{array}{cc}{S}_{{\mathrm{SSl}}_Y}\left(r_k,\stackrel{.fwdarw.}{r}\right)=\frac{\left(\stackrel{.fwdarw.}{r}-\stackrel{.fwdarw.}{.Math.\stackrel{.fwdarw.}{r}.Math.}\right)\stackrel{.fwdarw.}{.Math.l_Y.Math.}\left(r_k,\stackrel{.fwdarw.}{r}\right)}{\left(\stackrel{.fwdarw.}{\tau }-\stackrel{.fwdarw.}{.Math.\stackrel{.fwdarw.}{r}.Math.}\right)\stackrel{.fwdarw.}{.Math.L_Y.Math.}\left(\stackrel{.fwdarw.}{r}\right)}=\frac{\Delta {\mu }_{\mathrm{aY}}\left(\stackrel{.fwdarw.}{r}\right)}{\Delta {\mu }_{\mathrm{ak}}}& \left(16\right)\end{array}$

The sensitivity of the slopes of the three data types can be either positive or negative. This can be seen for the simple case of two source-detector separations, r.sub.1 and r.sub.2 (with r.sub.2>r.sub.1), and one focal change at r.sub.k. In this case, equation (16) becomes:

[00017]$\begin{array}{cc}{S}_{{\mathrm{SSl}}_Y}\left(r_k,\stackrel{.fwdarw.}{r}\right)=\frac{.Math.l_Y.Math.\left(r_k,r_2\right)-.Math.l_Y.Math.\left(r_k,r_1\right)}{.Math.L_Y.Math.\left(r_2\right)-.Math.L_Y.Math.\left(r_1\right)}& \left(17\right)\end{array}$

The interpretation of equation (17) is easier for the DC or AC data type. In this case, for focal changes of the absorption coefficient, the denominator is always positive while the numerator can be either positive or negative. For example, if a focal change occurs in proximity of the detector at r.sub.1 (where l.sub.Y(r.sub.k, r.sub.1)>l.sub.Y(r.sub.k, r.sub.2)), the sensitivity becomes negative. If a focal change occurs in proximity of r.sub.2 or deeper in the tissue, the sensitivity is positive. Also, there are deep regions of the medium where the slope method has a higher sensitivity to focal absorption changes than that of the corresponding data type (equation (10)). In fact, in this case l.sub.Y(r.sub.k, r.sub.2)>>l.sub.Y(r.sub.k, r.sub.i) and equation (17) has almost the same numerator of equation (10) but a smaller denominator. In general, the meaning of equation (16) and (17), at least for DC and AC data types, is clear: all the detectors located at distances less than (f)) give a negative contribution to the sensitivity, while all the detectors located at distances larger than (f)) give a positive contribution to the sensitivity. The reason for this property is that the partial and total pathlengths of these two data types are always positive, as discussed in § 2.1. Which contribution prevails, or whether they compensate each other, depends on the location and size of the region where a change in absorption occurs. For layered changes in the absorption coefficient, the two contributions compensate each other for a superficial layer location, yielding an almost null sensitivity. On the contrary, the contributions add up to be positive and bigger than the sensitivity of the corresponding data type for deeper location of the layer.

[0078] To correct, at least partially for the drawback of negative sensitivities in the single-slope method, the method described herein uses a special source-detector arrangement (typically, but not necessarily symmetrical). With this special source-detector array, which is described in § 2.4, it is possible to measure two matched single-slopes and take their average. This approach is referred to herein as a “dual-slope” method. The equations of this section are general, regardless of the geometry of the medium and heterogeneity of the optical properties. The only assumption made (beside the validity of first order perturbation of the diffusion equation) is that the functions y(r) are linear. In a practical situation, the equivalent absorption changes are found by inverting equation (14), which requires the knowledge of total pathlengths, as already discussed in the previous section.

2.3 Spatial Sensitivity to Focal Changes in Oxy- and Deoxyhemoglobin Concentrations

[0079] By reframing the changes of a data type (or its slope) to focal absorption perturbations with generalized pathlengths, it is possible to define the sensitivities to focal changes in chromophore concentrations. In the case of near-infrared spectroscopy of blood-perfused tissue, oxy- and deoxyhemoglobin are the dominant chromophores. Equations (9) and (15) can be written for two wavelengths (λ.sub.1 and λ.sub.2) to obtain two new equations that relate equivalent homogeneous changes in oxy- and deoxyhemoglobin (O and D, respectively) concentrations to their true focal changes. More precisely, the left sides of the new equations contain a linear combination of O.sub.Y and D.sub.Y (or O.sub.SSl.sub.Y and D.sub.SSl.sub.Y) i.e., the equivalent changes of oxy- and deoxyhemoglobin for the Y data type (or the slope of the Y data type), and the right sides contain a linear combination of the true focal hemoglobin changes (O.sub.k and D.sub.k). If S.sub.Y(λ.sub.1)≈S.sub.Y(λ.sub.2) (for the raw data Y) and S.sub.SSl.sub.Y(λ.sub.1)≈S.sub.SSl.sub.Y(λ.sub.2) (for the slopes of the data Y), equation (10) and equation (16) represent, with good approximation, also the sensitivity for a single focal change in hemoglobin concentration at r.sub.k (O.sub.k or D.sub.k) estimated with the data type Y or its slope, respectively. In other words, S.sub.Y≈O.sub.Y/O.sub.k≈D.sub.Y/D.sub.k, and S.sub.SSl.sub.Y≈O.sub.SSl.sub.Y/O.sub.k≈D.sub.SSl.sub.Y/D.sub.k The condition of the similarity of the absorption sensitivities at two wavelengths for equation (10) or equation (16) implies that the cross talk between the equivalent hemoglobin species (O.sub.Y and D.sub.Y, or O.sub.SSl.sub.Y and D.sub.SSl.sub.Y) is negligible.

2.4 Dual-Slope Method

[0080] A dual-slope method considers the average of two matched slopes obtained with a specially configured array of sources and detectors. Such a configuration features insensitivity to source and detector drift and to changes in opto-mechanical coupling between probe and tissue. These dual slopes offer some key advantages with respect to standard single-source or single-detector multi-distance methods used in near-infrared spectroscopy. In particular, dual slopes are less sensitive to localized superficial tissue inhomogeneities. Additionally, dual slopes feature a localized sensitivity to a deep tissue volume, as opposed to the typical banana-shaped region-of-sensitivity common in near-infrared spectroscopy.

[0081] The formulas presented herein are developed for a source-detector arrangement comprising first and second sources 20, 22 and several detectors 28-40 placed such that their distances from the two sources are pair-wise symmetric with respect to the average source-detector distance {right arrow over (r)} (FIG. 1D). The procedure used for deriving the formulas for the DC slope can also be used for AC and phase slopes. Application of the DC single-slope method (where y=ln[r.sup.2DC(r)] according to the previous section) for the first source 20 yields:

[00018]$\begin{array}{cc}{\mathrm{SSl}}_{\mathrm{DC}1}=\frac{\left(\stackrel{.fwdarw.}{r}-\stackrel{.fwdarw.}{.Math.\stackrel{.fwdarw.}{r}.Math.}\right)\stackrel{.fwdarw.}{\ln \left(r^2{\mathrm{DC}}_1\left(r\right)\right)}}{N\mathrm{var}\left(\stackrel{.fwdarw.}{r}\right)}=\frac{\underset{i=1}{\overset{N}{.Math.}}\left(r_t-.Math.\stackrel{.fwdarw.}{r}.Math.\right)\ln \left(r_i^2{\mathrm{DC}}_1\left(r_i\right)\right)}{N\mathrm{var}\left(\stackrel{.fwdarw.}{r}\right)}& \left(18\right)\end{array}$

where SSl.sub.DC1 is the DC single-slope calculated with respect to the first source 20 and N is the number of detectors 28-40. The formulas can be derived for either even N or odd N. However, for odd N, one of the detectors, such as the fourth detector 34 located at {right arrow over (r)}=20 millimeters in FIG. 1D, is placed at {right arrow over (r)}. This detector 34 does not contribute to the slope values or to its changes in time, therefore, it is useful to derive the formula for even N, in which case: r.sub.N-k−r=−[r.sub.k+1−r] with k=0, 1, . . . , N/2−1. By using the properties of logarithms, equation (18) becomes:

[00019]$\begin{array}{cc}{\mathrm{SSl}}_{\mathrm{DC}1}=\frac{\underset{k=0}{\overset{N/2-1}{.Math.}}\left(r_{N-k}-.Math.\stackrel{.fwdarw.}{r}.Math.\right)\left\{2\ln \left(\frac{r_{N-k}^{\prime }}{r_{k+1}^{\prime }}\right)+\ln \left(\frac{{\mathrm{DC}}_1\left(r_{N-k}\right)}{{\mathrm{DC}}_1\left(r_{k+1}\right)}\right)\right\}}{N\mathrm{var}\left(\stackrel{.fwdarw.}{r}\right)}& \left(19\right)\end{array}$

Similarly, for second source 22:

[00020]$\begin{array}{cc}{\mathrm{SSl}}_{\mathrm{DC}1}=\frac{\underset{k=0}{\overset{N/2-1}{.Math.}}\left(r_{N-k}^{\prime }-.Math.{\stackrel{.fwdarw.}{r}}^{\prime }.Math.\right)\left\{2\ln \left(\frac{r_{N-k}^{\prime }}{r_{k+1}^{\prime }}\right)+\ln \left(\frac{{\mathrm{DC}}_2\left(r_{N-k}^{\prime }\right)}{{\mathrm{DC}}_2\left(r_{k+1}^{\prime }\right)}\right)\right\}}{N\mathrm{var}\left({\stackrel{.fwdarw.}{r}}^{\prime }\right)}& \left(20\right)\end{array}$

[0082] In equation (20), the prime indicates distances calculated with respect to the second source 22. In FIG. 1D, it is assumed that r.sub.k=r′.sub.k which means for example, that the distance between the first source 20 and the second detector 30, is the same as the distance between the second source 22 and sixth detector 38. However, this is not strictly necessary for the validity of the dual-slope method, which can be applied under more general conditions. In a typical experiment, the measured DC signal from the first source 20 includes contributions related to source emission, detector response, and probe-to-tissue coupling:

[00021]$\begin{array}{cc}{\mathrm{DC}}_1\left(r_{N-k}\right)=P_1{C}_{1}Z\left(r_{N-k}\right)C_Z\left(r_{N-k}\right){\mathrm{dc}}_1\left(r_{N-k}\right)& \left(21\right)\\ {\mathrm{DC}}_1\left(r_{k+1}\right)=P_1{C}_{1}Z\left(r_{k+1}\right)C_Z\left(r_{k+1}\right){\mathrm{dc}}_1\left(r_{k+1}\right)& \left(22\right)\end{array}$

where P.sub.1 is the power emitted by the first source 20, C.sub.1 is the optical coupling factor between the first source 20 and tissue, Z(r.sub.N-k) [Z(r.sub.k+1)] is the sensitivity of the detector Z at r.sub.N-k [r.sub.k+1], and C.sub.Z(r.sub.N-k) [C.sub.Z(r.sub.k+1)] is the optical coupling factor between the detector Z at r.sub.N-k [r.sub.k+1] and tissue, and dc.sub.1(r.sub.N-k) and dc.sub.1(r.sub.k+1) are the theoretical values of DC intensities at a longer (r.sub.N-k) and shorter (r.sub.k+1) distance. The factors C.sub.1, Z, and C.sub.Z account for random or systematic temporal fluctuations, drifts in source or detector characteristics, displacement of the optical probe, etc.; whereas, dc.sub.1 does not include any kind of noise or experimental confounds. Similarly, for the second source 22:

[00022]$\begin{array}{cc}{\mathrm{DC}}_2\left(r_{N-k}^{\prime }\right)=P_2{C}_2{Z}^{\prime }\left(r_{N-k}^{\prime }\right)C_{Z^{\prime }}^{\prime }\left(r_{N-k}^{\prime }\right){\mathrm{dc}}_2\left(r_{N-k}^{\prime }\right)& \left(23\right)\\ {\mathrm{DC}}_2\left(r_{k+1}^{\prime }\right)=P_2{C}_2{Z}^{\prime }\left(r_{k+1}^{\prime }\right)C_{Z^{\prime }}^{\prime }\left(r_{k+1}^{\prime }\right){\mathrm{dc}}_2\left(r_{k+1}^{\prime }\right)& \left(24\right)\end{array}$

Because of the symmetrical source-detector arrangement, it is possible to rewrite equations (23) and (24) as:

[00023]$\begin{array}{cc}{\mathrm{DC}}_2\left(r_{N-k}^{\prime }\right)=P_2{C}_{2}Z\left(r_{k+1}\right)C_Z\left(r_{k+1}\right){\mathrm{dc}}_2\left(r_{N-k}^{\prime }\right)& \left(25\right)\\ {\mathrm{DC}}_2\left(r_{k+1}^{\prime }\right)=P_2{C}_{2}Z\left(r_{N-k}\right)C_Z\left(r_{N-k}\right){\mathrm{dc}}_2\left(r_{k+1}^{\prime }\right)& \left(26\right)\end{array}$

Averaging the single-slopes SSl.sub.DC1 and SSl.sub.DC2 and substituting equations (21) and (22) in equation (19) and equations (25) and (26) in equation (20) causes all the terms affected by temporal fluctuations to cancel out, thus leaving:

[00024]$\begin{array}{cc}\frac{{\mathrm{SSl}}_{\mathrm{DC}1}+{\mathrm{SSl}}_{\mathrm{DC}2}}{2}=\frac{\underset{k=0}{\overset{N/2-1}{.Math.}}\left(r_{N-k}-.Math.\stackrel{.fwdarw.}{r}.Math.\right)\left\{2\ln \left(\frac{r_{N-k}}{r_{k+1}}\right)+\frac{1}{2}\ln \left(\frac{{\mathrm{dc}}_1\left(r_{N-k}\right){\mathrm{dc}}_2\left(r_{N-k}^{\prime }\right)}{{\mathrm{dc}}_1\left(r_{k+1}\right){\mathrm{dc}}_2\left(r_{k+1}^{\prime }\right)}\right)\right\}}{N\mathrm{var}\left(\stackrel{.fwdarw.}{r}\right)}& \left(27\right)\end{array}$

where it has been assumed that r.sub.k=r′.sub.k, but no assumption has been made on the homogeneity of the medium.

[0083] Equation (27) defines an average “true” slope that depends only on the distribution of the optical properties in the medium (and the position of the source and the detectors) and that is not affected by any instrumental characteristics (laser power, detector sensitivity, coupling). For the case of a homogeneous semi-infinite medium, equation (27) defines a true slope that, in principle, could be calculated by using only a single source and N detectors under ideal conditions (i.e. same detectors sensitivities, same detectors coupling).

[0084] Under general conditions the dual-slope is defined as the average of the slopes obtained from the first source 20 and the second source 22:

[00025]$\begin{array}{cc}{\mathrm{DSl}}_{\mathrm{DC}}=\frac{{\mathrm{SSl}}_{\mathrm{DC}1}+{\mathrm{SSl}}_{\mathrm{DC}2}}{2}& \left(28\right)\end{array}$

Similarly, the dual-slopes of AC and phase are defined as:

[00026]$\begin{array}{cc}{\mathrm{DSl}}_{\mathrm{AC}}=\frac{{\mathrm{SSl}}_{\mathrm{AC}1}+{\mathrm{SSl}}_{\mathrm{AC}2}}{2}& \left(29\right)\\ {\mathrm{DSl}}_{\phi }=\frac{{\mathrm{SSl}}_{\mathrm{\phi 1}}+{\mathrm{SSl}}_{\mathrm{\phi 2}}}{2}& \left(30\right)\end{array}$

The sensitivity associated with the dual-slope method is:

[00027]$\begin{array}{cc}{S}_{{\mathrm{DSl}}_Y}=\frac{s_{{\mathrm{SSl}}_{Y1}+}{s}_{{\mathrm{SSl}}_{Y2}}}{2}& \left(31\right)\end{array}$

where S.sub.SSl.sub.Yi is the sensitivity associated with the single-slope calculated with respect to the source “i” (equation (16)).

[0085] Equation (31) implies also that the equivalent absorption change derived with the dual-slope method is the average of the equivalent absorption changes obtained with the two single-slopes (equation (15)). These source-detector arrangements not only achieve a signal that is not affected by source power, detector sensitivity, and optical coupling fluctuations, but also partially correct for the negative sensitivity of single-slopes present in some regions of the diffusive medium, as pointed out in § 2.2. In fact, if one considers a focal region close to the medium boundary and located in between the first source 20 and the first detector 24 (FIG. 1B), this region has a negative sensitivity with respect to the single-slope calculated by using the first source 20 and the first and second detectors 24, 26, but a positive sensitivity with respect to the symmetric single-slope, which uses the second source 22 and the same two detectors. By averaging the two single-slopes according to the dual-slope method, one achieves a partial compensation of the negative sensitivity. Finally, the dual-slope can be defined under more general conditions. In particular it is not required that r.sub.k=r′.sub.k. It is only required that that r.sub.k−{right arrow over (r)}=r′.sub.k−{right arrow over (r)}. This means that the second source 22 can be moved closer or further from the seventh detector 40 in FIG. 1D.

2.5 Solution of the Diffusion Equation in the Semi-Infinite Geometry

[0086] The foregoing equations define the sensitivity of the raw data at a single source-detector separation (equation (10)) and the slope of the raw data (equation (16)), the latter for both single-slope and dual-slope methods, for a semi-infinite medium geometry. For the single-distance data, results are shown only for the sensitivity at the furthest distance. In this geometry, the calculations of the partial pathlengths (equation (5)) are based on the following expression:

[00028]$\begin{array}{cc}\tilde{\overline{\Phi }}\left(r_s,r_k,\omega \right)=\frac{1}{4\pi D}\left[\frac{r_1}{}-\frac{r_2}{}\right]& \left(32\right)\end{array}$

which is the solution of the frequency-domain diffusion equation with extrapolated boundary condition for the fluence. The solution is calculated for a point-like source located at r.sub.s=(x.sub.0, 0,0) (where x.sub.0=1/μ′.sub.s; μ′.sub.s is the reduced scattering coefficient) and calculated at the point r.sub.k=(x.sub.k, y.sub.k, z.sub.k) inside the diffusive medium. The extrapolated boundary is located at x.sub.b=−2AD, where A is a parameter that considers the refractive index mismatch between diffusive medium and the outer medium, and D is the diffusion factor (D=1/(3μ′.sub.s)). By the method of images for solving partial differential equations, the virtual source is located at r′.sub.s=(−x.sub.0+2x.sub.b, 0,0), therefore r.sub.1=|r.sub.k−r.sub.s| and r.sub.2=|r.sub.k−r′.sub.s|. The complex effective coefficient is

[00029]$=\sqrt{\frac{{\mu }_a}{D}-i\frac{\omega }{\mathrm{vD}}}.$

[0087] Applying Fick's law to equation (32) yields the complex reflectance (the output complex intensity at r.sub.d when photons are emitted from r.sub.s):

[00030]$\begin{array}{cc}\tilde{R}\left(r_{s\prime }{r}_d\omega \right)=\frac{1}{4\pi }\left[x_0\left(\frac{1}{r_1}+\right)\frac{r_1^2}{}+\left(x_0-2x_b\right)\left(\frac{1}{r_2}+\right)\frac{r_2^2}{}\right]& \left(33\right)\end{array}$

where, in this case, r.sub.1=|r.sub.d−r.sub.s|, r.sub.2=|r.sub.d−r′.sub.s|. As for {tilde over (R)}(r.sub.k, r.sub.d, ω), equation (33) is used with x.sub.k replacing x.sub.0 and r.sub.1=|r.sub.d−r.sub.k|, r.sub.2=|r.sub.d−r′.sub.k|. In fact, for this case the real and virtual sources are r.sub.k and r′.sub.k=(−x.sub.k+2x.sub.b, y.sub.k, z.sub.k), respectively. The complex total pathlength {tilde over (L)}(r.sub.d, ω) is calculated as:

[00031]$\begin{array}{cc}.Math.\tilde{L}.Math.\left(r_s,r_{d\prime }\omega \right)=\frac{1}{8\pi D\tilde{R}\left(r_s,r_d,\omega \right)}\left[\left(\frac{x_0}{r_1}\right)+\left(\frac{x_0-2x_b}{r_2}\right)\right]& \left(34\right)\end{array}$

where r.sub.1 and r.sub.2 have the same meaning as in equation (33). Equations (32)-(34) are valid for all the data types in the frequency-domain and reduce to the DC expressions when ω=0. The calculations of the mean partial pathlengths were carried out by numerical integration, by dividing a region into smaller voxels of 1 mm.sup.3 volume. All frequency-domain results were obtained at a modulation frequency f=140 MHz (ω=2πf). The geometrical arrangement of sources and detectors is illustrated in FIG. 1C. A particularly suitable frequency of modulation at about 100 MHz.

2.6 Considerations on Noise-to-Signal Ratio

[0088] Thus far, theoretical results have been derived without considering the noise affecting different data types and the associated slopes. This section takes noise into account.

[0089] As used herein, ΔY.sub.noise refers to the noise-to-signal ratio for the data type Y. For DC data,

[00032]$\Delta {\overline{Y}}_{noise}=\frac{{\sigma }_{DC}}{DC}.$

where σ.sub.DC is the standard deviation and DC is the average of DC data calculated in a given time interval (similar formula applies for AC). For phase data ΔY.sub.noise=σ.sub.Ph, where σ.sub.Ph is the standard deviation of phase. Typical values of

[00033]$\frac{{\sigma }_{DC}}{DC},\frac{{\sigma }_{AC}}{AC},$

and σ.sub.Ph are 0.1%, 0.1%, and 0.1°, respectively. For the single-slope where only two distances are used (r.sub.1, r.sub.2; r.sub.2>r.sub.1), the changes in slope are given by:

[00034]$\begin{array}{cc}\Delta {\mathrm{SSl}}_Y=\frac{\Delta \overline{Y}\left(r_2\right)-\Delta \overline{Y}\left(r_1\right)}{r_2-r_1}& \left(35\right)\end{array}$

By the formula of a priori error propagation, the following equation gives the noise level of a slope change:

[00035]$\begin{array}{cc}\Delta SS{l}_{\mathrm{Y\_noise}}\times \left(r_2-r_1\right)=2\Delta {\overline{Y}}_{noise}& \left(36\right)\end{array}$

where ΔSSl.sub.Y_noise indicates the noise level of slope change for data type Y and where there exists equal noise level at the two source-detector separations. A better representation of the error in the slopes is obtained by adding the errors of the raw data at the two distances quadratically to yield (for the case of equal noise-to-signal ratio at both distances):

[00036]$\begin{array}{cc}\Delta SS{l}_{\mathrm{Y\_noise}}\times \left(r_2-r_1\right)=\sqrt{2}\Delta {\overline{Y}}_{noise}& \left(37\right)\end{array}$

[0090] Equation (37) assumes that the raw data does not contain relevant drifts. If it does, error propagation will yield an overestimation of the errors in the slopes. The same is true for the error in the dual-slope method. Adding the errors in the single-slopes quadratically yields the correct representation of the error in the dual-slope only if any relevant drift is absent in the two slopes. Making this assumption and assuming equal errors in the two slopes results in:

[00037]$\begin{array}{cc}\Delta DS{l}_{\mathrm{Y\_noise}}\times \left(r_2-r_1\right)=\Delta {\overline{Y}}_{noise}& \left(38\right)\end{array}$

where ΔDSl.sub.Y_noise is the noise level of dual-slope change for data type Y. Using the relationship between change in a data type (due to a focal change in absorption at r.sub.k) and the corresponding partial pathlengths results in a useful way to determine a limit of detectability in simulated data (i.e. when the source of contrast is equal to the noise level). For the single-distance data Y:

[00038]$\begin{array}{cc}\Delta {\overline{Y}}_{noise}=\left|.Math.l_Y.Math.\left(r_k,r\right)\Delta {\mu }_a\right|& \left(39\right)\end{array}$

where the left side is equal to 0.1% for DC or AC data and 1.7.Math.10.sup.−3 rad(0.1°) for phase data.

[0091] For the single-slope method:

[00039]$\begin{array}{cc}\Delta SS{l}_{\mathrm{Y\_noise}}\times \left(r_2-r_1\right)=\left|\left[.Math.l_Y.Math.\left(r_k,r_2\right)-.Math.l_Y.Math.\left(r_k,r_1\right)\right]\Delta {\mu }_a\right|& \left(40\right)\end{array}$

where l.sub.Y(r.sub.k, r.sub.i) is the pathlength spent by detected photons (at the source-detector separation r.sub.i) inside the focal region at r.sub.k. The left side is equal to √{square root over (2)}.Math.0.1% for DC and AC data and √{square root over (2)}.Math.1.7.Math.10.sup.−3 rad for phase data. For the dual-slope method:

[00040]$\begin{array}{cc}\Delta DS{l}_{\mathrm{Y\_noise}}\times \left(r_2-r_1\right)=\frac{1}{2}\left|\left\{\left[{.Math.l_Y.Math.}_{S_1}\left(r_k,r_2\right)+{.Math.l_Y.Math.}_{S_2}\left(r_k,r_2\right)\right]-\left[{.Math.l_Y.Math.}_{S_1}\left(r_k,r_1\right)+{.Math.l_Y.Math.}_{S_2}\left(r_{k\prime }{r}_1\right)\right]\right\}\Delta {\mu }_a\right|& \left(41\right)\end{array}$

where l.sub.Ys.sub.j(r.sub.k, r.sub.i) indicates the partial pathlength (for the data type Y) inside the focal defect at r.sub.k when photons are detected at a source-detector separation r.sub.i calculated from the source 20, 22. The left-hand side is equal to 0.1% for DC or AC data and ˜0.0017 rad(0.1°) for phase data. In equations (39)-(41) the right-hand sides are known from the simulations and they are the absolute values of the source of contrast. The limits of detectability are defined when the source of contrast is equal to the noise level (left side of equations (39)-(41)). By using equations (10) and (16) it is possible to also determine the noise level of the equivalent absorption change and the noise level of oxy- and deoxyhemoglobin oscillations. The noise level of equivalent absorption change is wavelength dependent, therefore in a given scenario where oxy- and deoxyhemoglobin oscillations are associated with different tissue regions, it is possible that only the equivalent absorption change at one wavelength is above noise level. This may happen in situations when the noise-to-signal ratio is approximately unity. In this case nothing can be concluded about the equivalent oxy- and deoxyhemoglobin changes derived with the method proposed (i.e. if they are below or above noise level), unless the errors are propagated to determine their noise levels. In the following discussions equations (39)-(41) are applied to determine the noise level for the change of a data type (or the change of its slope) in numerical results.

3. Results

[0092] The results are presented in four sections, one for layered absorption changes (§ 3.1), two for focal absorption changes (Sections 3.2 and 3.3), and one where that included application of noise-to-signal ratio considerations to numerical results (§ 3.4). Both layered and focal changes are of interest for describing hemodynamics that are usually studied in functional near-infrared spectroscopy of the human head. When one studies cerebral hemodynamics of systemic origin, as is the case in coherent hemodynamics spectroscopy, layered absorption changes are often a good approximation.

[0093] For typical cases of brain activation, the hemodynamic changes are usually more complex. In such cases, focal changes occurring in the brain cortex and hemodynamic activity that occurs in the extracerebral tissue layers tends to confound measurements. In each case, the sensitivity of single-distance data (equation (10)), single-slope data (equation (16)), and dual-slope data (equation (31)). For the single-distance data, the farthest distance is assumed (35 millimeters).

[0094] There exist different methods to recover periodic oscillations in oxy- and deoxyhemoglobin, which are described by the phasors O and D, respectively. Such methods are described in § 3.1 for layered changes and in § 3.3 for focal changes, respectively.

[0095] In particular, in coherent hemodynamics spectroscopy, it is of interest to measure D/O, i.e. the ratio of amplitudes and the phase difference between the two hemoglobin species. For the simulations, two (true) phasors for O and D (O.sub.T and D.sub.T) are assigned to each region of the medium where oscillations occurred. The phasors O.sub.T and D.sub.T are then translated into (true) absorption phasors at two wavelengths (690 and 830 nanometers). The phasors are then recovered corresponding to the equivalent absorption change for the single-distance data, for the single-slope and dual-slope by using equation (9), equation (15), and the equation derived from equation (31), respectively. The recovered phasors of absorption oscillations are translated into estimated phasors for oxy- and deoxyhemoglobin O.sub.E and D.sub.E. Finally, D.sub.T/O.sub.T (defined for each region) is compared with D.sub.E/O.sub.E, which is the equivalent phasor ratio obtained for each method. § 3.4 revisits the numerical results of the previous three sections by adding a few comments on noise-to-signal ratio.

3.1 Sensitivity to Layered Absorption Changes

[0096] One example of sensitivity for layered changes in the absorption coefficient is shown in FIG. 2, obtained with the source-detector arrangement of FIG. 1C. The sensitivities were obtained according to equation (10) (FIG. 2, top panel) or equation (16) (FIG. 2, bottom panel) by scanning a layer of size 1×80×80 millimeters (in the x, y, and z directions), along the depth of the medium (x axis) with 1 millimeter increments. The center of the layer was fixed in the (y and z) directions at (20 and 0) millimeters. The optical properties of the medium were: μ.sub.a=0.01 mm.sup.−1 and μ′.sub.s=1 mm.sup.−1; the refractive indices of the diffusive medium and outer medium were: n.sub.i=1.4 and n.sub.o=1, respectively. The modulation frequency for AC and phase data was 140 MHz. For this case of layered changes, the results of the slope method are the same regardless of whether the single-slope is being considered (with respect to the first source 20 or the second source 22 of FIG. 1C) whether they are averaged according to the dual-slope approach. Therefore in this case both sensitivities (with single- or dual-slope) are indicated as S.sub.Sl.sub.Y. Two important features of FIG. 2 are: (1) the shift towards the right (deeper regions of the medium) of the maxima of S.sub.Sl.sub.Y with respect to the maxima of S.sub.Y; (2) the almost null sensitivity of S.sub.Sl.sub.Y when the layer is close to the medium boundary. This result is because the region of the layer close to the first detector 24 (close detector) has a negative sensitivity which is compensated by the positive sensitivity of the region close to the second detector 26 (far detector). The reason why DC and AC do not have the maximum sensitivity when the layer is at x=0.5 millimeters (top position), is because these data types have a high sensitivity mostly around source and detector, but very little in the central region of the layer (this region is crossed mainly by few photons with short paths). As one scans the layer further down to a depth of 5 millimeters, the sensitivity to the layer region close to source and detector is smaller, but it is bigger in the central part of the layer (now more photons with longer paths will cross the entire layer), causing an overall increase in sensitivity.

[0097] The negligible value of S.sub.Sl.sub.Y to absorption changes occurring at shallow depths shows that a slope method offers a natural “depth filter” for cases when hemodynamics oscillations occurs both close to the boundary and deeper in the medium. We prove this point by considering two cases of DSO estimation. The first case considers hemodynamic oscillations that occur only in two layers at the same time: a top layer 1 millimeter thick which occupies the region xϵ[0,1] millimeters, and a deeper layer 1 millimeter thick which occupies the region xϵ[15,16] millimeters. These layers are representative of a skin layer (or top layer “T”) and a brain layer (or bottom layer “B”). The phasors describing the true oscillations are: O.sub.1=O.sub.2=1e.sup.i0 μM, D.sub.1=0.4e.sup.i0 μM, and D.sub.2=0.4e.sup.i3/4π μM. The subscripts “1” and “2” indicate the top layer and bottom layer, respectively. Therefore:

[00041]$\frac{D_{1T}}{O_{1T}}=0.4e^{i0}\mathrm{and}\frac{D_{2T}}{O_{2T}}=0.4e^{i3/4\pi }.$

[0098] FIG. 3 compares the estimated ratio of phasors D.sub.E/O.sub.E with DC and DC slope method, which is shown in the left panel, with that obtained using phase and phase slope method, which is shown in the right panel. The labels “T” and “B” indicate the true ratio of phasors relative to top and bottom layer D.sub.T/O.sub.T. The ratio of phasors estimated with AC and AC slope is almost identical to those obtained with DC and DC slope, respectively.

[0099] The left panel shows that while DC data yields a ratio of phasors “closer” to the true one of the top layer, the DC slope method yields a ratio of phasors very close to the true one in the bottom layer. Phase and the phase slope method yield phasor ratios that almost overlap with the true phasor ratio of the bottom layer.

[0100] The results shown in FIG. 3 were obtained using a modulation frequency of f=140 MHz and background optical properties of μ′.sub.s (690 nanometers)=1.2 mm.sup.−1 and μ′.sub.s (830 nanometers)=1 mm.sup.−1. The absorption coefficients at the two wavelengths were calculated by considering a total hemoglobin concentration ([HbT]) of 45 μM and a hemoglobin saturation (St) of 0.65. Similar results (not shown) were obtained for μ′.sub.s (690 nanometers)=1.5 mm.sup.−1 and μ′.sub.s (830 nanometers) 1.2 mm.sup.−1 and for the absorption coefficients at the two wavelengths relative to [HbT] of 75 μM and St of 0.65.

[0101] The second case considers hemodynamic oscillations that occur in two layers at the same time: a top layer 6 millimeters thick that occupies the region xϵ[0,6] millimeters, (indicated by “T”) and a deeper layer 6 millimeters thick that occupies the region xϵ[12,18] millimeters (indicated by “B”). The phasors describing the true oscillations (red arrows) are the same as before for the top and bottom layers. The background optical properties and modulation frequency are also the same as those of the previous case. The results are shown in FIG. 4.

[0102] The case shown in FIG. 4 may be a closer representation of what happens during induced hemodynamic oscillations in coherent hemodynamics spectroscopy, with systemic oxy- and deoxyhemoglobin changes occurring both in the scalp and in the brain, and most likely a region of negligible oscillations in the skull and cerebrospinal fluid (CSF) layers. As before, the ratios of phasors estimated with DC and DC slope are on the left panel, while those estimated with phase and phase slope are on the right panel. For this case, except the phasor's ratio estimated with phase slope, which is coincident with the true phasor's ratio of the bottom layer, all the other data types and slopes show a different degree of mixture between the true phasor ratio of top and bottom layer. The phasor ratio estimated with DC is the closest to the true phasor ratio of the top layer (left panel); this is followed by the one obtained with DC slope which points at an angle in between the two true phasor ratios (left panel). The phasor ratio estimated by phase is closer to the true one of the bottom layer (right panel). Changing the geometry of the two layers or the magnitude of oscillations (by keeping the true phasor ratio fixed), which in the real tissue might be stronger in the brain, results in all the phasor ratios estimated by the raw data and the slopes method would rotate counter clock wise toward the true phasor ratio of the bottom layer. It would have been possible to obtain similar results by directly applying the sensitivity functions equation (10) and (16) to the ratio of phasors describing oxy- and deoxyhemoglobin oscillations (O.sub.E/O.sub.T and D.sub.E/D.sub.T). For the case of FIG. 3, the discrepancy between the phasor ratios retrieved by using the two methods is less than 5° and less than 9% for the angles and the amplitudes, respectively.

3.2 Sensitivity to Focal Absorption Changes

[0103] This section concerns the study of the sensitivity of the single-distance data and their slopes to focal changes in absorption. Only DC and phase data is reported herein. This is because of the similarity of AC data and derived slopes with those of DC data. By considering the source-detector arrangement of FIG. 1C, it is possible to scan a rectangular cuboid of size 1×5×5 millimeters (in the x, y, and z directions, respectively) was scab along the y (horizontal) and x (vertical) directions at 1-millimeter step increments. The sensitivity maps reported in FIG. 5 were obtained by considering the source-detector pair that consists of the first source 20 and the second detector 26 in FIG. 1C (distance of 35 millimeters). As expected from § 2.1, DC always shows a positive sensitivity, while the phase has a more complex behavior. The computations were performed for the optical properties and modulation frequency are the same as those of FIGS. 3 and 4 at 690 nanometers. The white (saturated) pixels in the DC map represent values higher than approximately 4.5×10.sup.−3.

[0104] Both DC intensity and phase sensitivity maps have a “banana” shape. FIG. 6 shows cross-sectional lines of the map of FIG. 5 at different depths (x=0.5, 3.5, 6.5, 9.5, 12.5, 15.5 millimeters). Some notable features of FIG. 6 are: (a) DC always shows a positive sensitivity, while the phase can have either a positive or negative sensitivity (as described in § 2.1); (b) with respect to DC, phase data is less sensitive close to the medium boundary and more sensitive in deeper regions of the diffusive medium; (c) the asymmetric behavior of DC intensity at a depth of x=0.5 millimeters results from the introduction of an effective source inside the medium, breaking the symmetry of source-detector configuration. Also, the voxel discretization used for the calculation in this work caused more photons to be “absorbed” when the region was close to the detector than when it was close to the source. The asymmetry is not visible in phase data, because for focal absorption increases in both locations (i.e. close to the source and close to the detector) comparable fractions of short-path and long-path photons were “absorbed”, and the phase was affected almost equally. FIG. 7 shows the sensitivity maps for single-slopes of DC (top panel), and phase (bottom panel). The two sets of saturated pixels represent values less than −0.01 (black) and higher than 0.02 (white). They were obtained by using a first source-detector pair, which consists of the first source 20 and the first detector 24 and a second source-detector pair, which consists of the first source 20 and the second detector 26, as shown in FIG. 1C, with distances of 20 and 35 millimeters respectively. As discussed in § 2.2, S.sub.SSl.sub.Y can be either positive or negative. The locations of a focal region where S.sub.SSl.sub.Y becomes negative are in proximity of the detector at the shortest distance from the source (the first detector 24, which is 20 millimeters from the first source 20). Even though the slope method features a greater sensitivity to deeper tissue regions with respect to single-distance data, the areas of negative sensitivity can be problematic for the interpretation of the data. FIG. 8 represents some cross-sectional lines of the sensitivity maps of FIG. 7 for single-slopes, where each plot refers to a different depth. As shown in FIG. 8, the phase slope has a smaller sensitivity than DC slope at shallow depths (x=0.5 and 3.5 millimeters), but a greater sensitivity in deeper regions of the tissue (x>12.5 millimeters). All the single-slopes show positive as well as negative sensitivities, especially for absorption changes close to the first detector 24. For DC, the single-slope sensitivity is very large and negative in proximity of the first detector 24 (close detector) and very large and positive (even greater than for single-distance data) close to the second detector 26 (see plot at x=0.5 millimeters). By comparing the cross-sectional lines in FIGS. 6 and 8, it is apparent that deeper in the tissue (for example x=12.5 and 15.5 millimeters) the sensitivity of DC, and phase single-slopes are greater than for the corresponding single-distance data. The dual-slope method compensates, at least partially, for the drawbacks of the single-slopes, i.e., the negative sensitivity and the strong positive sensitivity to focal perturbations close to the detectors.

[0105] The dual-slope method uses the short distance between the first source 20 and the first detector 24 and the long distance between the first source 20 and the second detector 26 for calculating a first slope. It then uses the short distance between the second source 22 and the second detector 26 and the long distance between the second source 22 and the first detector 24 for calculating the second slope.

[0106] However, the result would not change if one were to use the short distance between the first source 20 and the first detector 24 and the long distance between the second source 22 and the first detector 24 for one slope and the short distance between the second source 22 and the second detector 26 and long distance between the first source 20 and the second detector 26 for the other slope.

[0107] The two slopes are averaged according to the dual-slope method. The sensitivity maps of DC and phase with the dual-slope method are shown in FIG. 9, in which the white pixels indicate regions where the sensitivity is higher than approximately 7×10.sup.−3. A comparison of FIGS. 9 and 7, confirms that the regions of negative sensitivity of the dual-slopes are smaller in size and feature a reduced sensitivity value with respect to the single-slopes. This feature is more obvious when one looks at the cross-sections of the sensitivity along the y axis at different depths, as shown in FIG. 10 (DC dual-slope) and FIG. 11 (phase dual-slope). From FIG. 9, one can see that the phase dual-slope has more desirable features than the DC intensity dual-slope, in that its maximal sensitivity to focal absorption changes is deeper in the tissue (about eleven millimeters as opposed to about five millimeters), and its sensitivity to superficial (about a two millimeter depth) absorption changes is lower (−0.004 compared to −0.0055).

[0108] A quantitative comparison between the sensitivities obtained with the dual-slope method and those obtained with single-distance data shows that the advantages of the former outweighs its only drawback, which is the presence of some residual negative values close to the medium boundary. When one compares DC dual-slope and single-distance data, one observes that: (a) the dual-slope DC has a much more spatially confined region of positive sensitivity than the typical “banana shape” of single-distance DC data (FIG. 5); (b) the DC dual-slope has a higher positive sensitivity than single-distance DC data (except for superficial tissue: x=0.5 millimeters). Similarly, if one compares the sensitivities obtained with dual phase slope and single-distance phase data one sees that: (a) the dual phase slope has a much more spatially confined region of positive sensitivity compared with the positive “banana” of the single-distance phase data (FIG. 5); (b) the dual phase slope has a smaller positive sensitivity than the single-distance phase data down to a depth of x=6.5 millimeters, and a higher sensitivity deeper in the medium; (c) the region of negative sensitivity is more spatially confined than that of the single-distance phase data, but the values are slightly more negative.

[0109] It is possible to further reduce the negative sensitivities present in the dual-slope method if one considers a multi-detector arrangement (FIG. 1D), as shown in FIG. 12. This arrangement is currently only possible with DC instrumentation that feature a suitably large dynamic range. Of course, one can choose other arrangements which use more than two detectors but require a smaller dynamic range. A comparison of FIGS. 9 and 12 shows that even though the negative sensitivities are smaller for the configuration featuring seven detectors, the region of positive sensitivity are also slightly smaller and less spatially confined than those obtained with two detectors.

3.3 Recovery of Localized Oscillations of Oxy- and Deoxyhemoglobin Concentrations

[0110] The following discussion considers two examples of phasors ratio retrieval (D/O) as it was carried out for the case of layered absorption changes. These examples reflect situations of periodic brain activation (in response to a block paradigm) where a focal change in cerebral hemodynamics occurs synchronously with hemodynamic changes in the extracerebral tissue layers. The hemodynamic changes in the brain are often larger than in the extracerebral layers.

[0111] The geometry is shown in FIG. 13. In both cases, what follows is a comparison of only the results obtained with the single-distance data and with the dual-slope methods because the single-slope methods yielded results that are almost coincident with those of the dual-slope methods. This is true because of the particular focal perturbations used (which involve absorption changes over a superficial layer with a relatively large area in the y-z plane).

[0112] In a first case, there are three top regions (6-millimeter thick and 22×40 millimeters in lateral dimensions) that are characterized by the same oscillations of oxy- and deoxyhemoglobin concentrations, expressed by the phasors O.sub.1=O.sub.2=O.sub.3=1e.sup.i0 μM, and D.sub.1=D.sub.2=D.sub.3=0.4e.sup.i0 μM, which represent a blood volume oscillation. It is also possible to consider a deeper region (a cube of 10-millimeter sides) that is characterized by the phasors O.sub.4=5e.sup.i0 μM, and D.sub.4=5e.sup.i3/4π μM, which represent an almost pure blood flow oscillation (for a pure blood flow oscillation deoxy- and oxyhemoglobin would shifted by it).

[0113] FIG. 14 shows the retrieval of the phasor ratio D.sub.E/O.sub.E by using single-distance DC and dual-slope DC (left panel), and single-distance phase and dual-slope phase (right panel). A second case (FIG. 15) considers different hemodynamics in the top three regions O.sub.1=O.sub.2=O.sub.3=1e.sup.i0 μM, D.sub.1=0.4e.sup.i0 μM, D.sub.2=0.4e.sup.iπ/3 μM, and D.sub.3=0.4e.sup.iπ/2 μM and in the deeper region O.sub.4=5e.sup.i0 μM and D.sub.4=5e.sup.i3/4π μM. As can be seen in FIGS. 14 and 15, the closest measure of the phasor ratio for the deeper region is obtained the dual phase slope method. The DC dual-slope and single-distance phase gave similar results. In both cases, the single-distance DC data always yielded a phasor ratio closer to one of the three top regions. Note that the amplitudes of the retrieved phasors with the different methods, O.sub.E, D.sub.E, are very similar to values measured in a typical coherent hemodynamics spectroscopy protocol.

3.4 Analysis of Signal to Noise Ratio: Application to Numerical Results

[0114] This section uses the equations describing the noise-to-signal ratio (equation (39)-(41)) to revisit the example of FIG. 15 and to add further considerations about the signal-to-noise ratio (SNR (λ), which is the inverse of noise-to-signal ratio). The purpose is twofold: (1) to learn why a retrieved ratio of phasors is in the direction it is found (2) to learn which region (characterized by certain phasors) in the medium is detectable (i.e. has SNR>1). This is done for a typical case of brain activation with some superficial confound (FIG. 15).

[0115] To answer the first question, one does not need to use the concept of SNR. It suffices, for each region, to calculate the source of contrast (right side of equations (39)-(41), i.e. the signal) and compare them. The source of contrast (or signal) in one region depends on the partial pathlength spent by detected photons in that region and by the change in the absorption coefficient. Since the oscillations of oxy- and deoxyhemoglobin in each region have been specified, it is possible to determine the contrast in each a region. If the source of contrast in one region is much bigger than that in the others, the retrieved ratio of phasors will be close to the true ratio of phasors in that region. If the sources of contrast in different regions are similar, the retrieved ratio of phasors will be pointing in an intermediate direction.

[0116] These considerations apply also if none of the regions are detectable. On the contrary, dividing the source of contrast by noise adds the information if a region is detectable or not. Several regions may not be detectable separately but may prove to be detectable when they are combined.

[0117] A final example reports the signal-to-noise ratio maps for changes of absorption coefficient derived from localized changes of oxy- and deoxyhemoglobin concentration where both the volume of the region and the changes of hemoglobin concentrations are typical of brain activation. The purpose is to learn which method can be used to detect a focal perturbation in a realistic case of brain activation.

[0118] SNR is defined as: SNR.sub.Ŷ(λ)=ΔY(λ)/ΔY.sub.noise (λ) for the single-distance data Y, SNR.sub.SSl.sub.Y(λ)=ΔSSl.sub.Y(λ)/ΔSSl.sub.Y_noise(λ) for the single-slope, and SNR.sub.DSl.sub.Y(λ)=ΔDSl.sub.Y(λ)/ΔDSl.sub.Y_noise(λ) for the dual-slope. More precisely ΔY, ΔSSl.sub.Y, and ΔDSl.sub.Y (the signal levels) are given by the right-hand sides of equations (39)-(41), respectively, and ΔY.sub.noise, ΔSSl.sub.Y_noise, and ΔDSl.sub.Y_noise (the noise levels) are those discussed in § 2.6. The simulations use only two wavelengths: λ.sub.1=690 nanometers and λ.sub.2=830 nanometers. For each data type Y, the SNR analysis is carried out only at the farthest source-detector separation (and specifically for the first source 20 and the second detector 26 of FIG. 13). For the single-slopes, only source the first source 20 and first and second detectors 24, 26 A and B of FIG. 13 were considered. In the example of FIG. 15, different oscillations of oxy- and deoxyhemoglobin concentrations were assigned to different focal regions of the medium and assumed to be present at the same time. Even in this case, the proposed method always retrieves one equivalent oscillation for oxy- and deoxyhemoglobin (for each single-distance data type or associated single-slope or dual-slope case), which is the results of the combined effect of all the concentration oscillations in the medium. However, by calculating the right-hand sides of equations (39) and (41) (which represent the signal level) separately in each region where oscillations occur, it is possible to calculate their SNR. This provides information on how closely the retrieved equivalent oscillations will reproduce the real ones in each region.

[0119] The SNR results for the example of FIG. 15 are shown in FIG. 18 for single-distance and dual-slope data (for both DC intensity and phase) associated with absorption changes in regions 1, 2, 3, and 4 (4 being the cubic region). With single-distance and dual-slope phase data, only the cubic region 4 yielded a SNR above the noise level. With single-distance phase, SNR.sub.φ(λ.sub.1)=3.5 and SNR.sub.φ(λ.sub.2)=1.3; however, these two SNR.sub.φ accounts for 77% and 56% of the total SNR.sub.φ at λ.sub.1 and λ.sub.2, respectively. In fact, even if no other region yielded a SNR above noise when considered individually, the combined effect of regions 1, 2, 3 did. For these reasons the retrieved phasors ratio with single-distance phase is not exactly coincident with the true phasor ratio of the cubic region. With dual phase slope SNR.sub.DSl.sub.φ(λ.sub.1)=3.3, SNR.sub.DSl.sub.φ(λ.sub.2)=1.2, which account for 95% and 91% of the total SNR, respectively. Therefore, the phasor ratio obtained with dual phase slope should be the closest to the true one of the cubic region, as it is shown in FIG. 15. These numbers give some indication about the retrieved phasor ratio with the different methods. For example, for DC single-distance data, regions 1, 2, and 4 contribute almost equally to the total SNR at λ.sub.1(33%, 32%, and 35% respectively), while at λ.sub.2 the SNR originates mostly from regions 1 and 2 (41% and 40%, respectively) and less from region 4 (18%). For the DC dual-slope, region 1 has a small contribution to the total SNR (4% at λ.sub.1). Slightly bigger contributions with respect to DC data comes from region 2 (34% and 54% at λ.sub.1 and λ.sub.2, respectively). The largest increase of SNR with respect to DC single-distance data achieved with DC dual-slope is found in the cubic region (59% and 37% at λ.sub.1 and λ.sub.2, respectively). These considerations are consistent with a ratio of phasor retrieved with DC dual-slope which is closer to the true ratio of phasors of the cubic region than that obtained with DC single-distance data (FIG. 15). Similar reasoning applies also for phase single-distance and phase dual-slope data (FIG. 15), as discussed above.

[0120] Of particular interest is the SNR for the case of a “realistic” brain activation. During visual stimulation, a local increase of about 30% in cerebral blood volume has been measured. We model a focal brain activation as a cuboid of sides 2×10×10 millimeters in the x, y, and z directions, where 2 millimeters is representative of the thickness of the brain cortex. The cuboid was scanned in a semi-infinite medium geometry with 1-millimeter step in both x and y directions. The background optical properties of FIG. 3) (at λ.sub.1=690 nanometers) were used. For the perturbation, [HbO] was increased by 9 μM and [Hb] was increased by 5 μM (which caused Δμ.sub.a(λ.sub.1)=0.0033 mm.sup.−1). This kept hemoglobin saturation (St) fixed at 0.65.

[0121] For consistency with the sensitivity maps (FIGS. 5, 7, and 9), equations (39)-(41) were used with the right terms without absolute value and with a minus sign preceding them. This defines a SNR with a sign with the properties that its absolute value represents the true SNR (positive by definition) and its sign is the same of the equivalent absorption change. The maps of SNR are reported for DC, DC single-slope and DC dual-slope (FIG. 16), and for phase, phase single-slope and phase dual-slope (FIG. 17).

[0122] FIGS. 16-17 highlight the following points: (1) in general DC data (single-distance, single-slope, dual-slope) are characterized by a much higher signal-to-noise ratio at all depths than phase data; (2) in particular, while DC single-distance data shows the highest signal-to-noise ratio at all depths, it is the highest in the top 4-5 millimeters and especially in regions close to the source and the detector (saturated white pixels in the top panel of FIG. 16), which makes this data type particularly sensitive to superficial changes in absorption; (3) DC single-slope has a smaller and more spatially localized SNR compared to single-distance DC data, but it also features a region of negative SNR (reflecting the negative sensitivity) close to the detector at the shorter distance (see middle panel of FIG. 16); (4) the signal obtained with DC dual-slope has the most spatially confined positive SNR (see bottom panel of FIG. 16) and the regions with negative SNR have smaller values when compared with single-slope. This is due to the partial cancellation of the negative sensitivities intrinsic in the dual-slope method; (5) phase data shows a “banana shaped” region of positive SNR (above noise level) with rather constant values down to depths of 12-13 millimeters. There is also a central top region (between the first source 20 and farthest detector) of negative SNR which is below the noise level, showing that phase data is not sensitive to absorption changes occurring close to the medium boundaries; (6) phase dual-slope has the most desirable features for imaging. In fact, not only does it show a negative signal-to-noise ratio below noise in the top 6-7 millimeters, it also shows a very confined region of positive SNR that spans the depths ˜8-15 millimeters and about 10-15 millimeters in width. The only drawback is the low signal-to-noise ratio.

[0123] Similar results were found also for different background properties characterized by [HbT] of 75 μM, St of 0.65, μ′.sub.s(λ.sub.1)=1.5 mm.sup.−1, and μ′.sub.s(λ.sub.2)=1.2 mm.sup.−1. In this case the 30% increase in [HbT] (by keeping a constant St) was obtained by increasing the oxy- and deoxyhemoglobin concentrations at the defect by 15 μM and 8 μM, respectively, relative to background values.

[0124] The preceding discussion shows that, through diffusion theory, dual slopes based on DC intensity, AC amplitude, or phase data types achieve a greater relative sensitivity to deep vs. shallow regions in a diffusive medium with respect to single slopes.

[0125] Compared to single-distance and single-slope data, the dual-slope method also features a more spatially confined region-of-sensitivity, which does not extend to the medium surface, as it is the case for banana-shaped regions-of-sensitivity, and reaches more deeply into the medium when phase rather than intensity data is used.

[0126] The methods and systems described herein rely at least in part on the concept of equivalent absorption change, which is defined as the homogeneous change in the absorption coefficient of the whole medium that yields the same data type change (or change in the slope of a data type) that is measured when one or more focal absorption changes are present in the medium. The definition of equivalent absorption change for each data type (and derived single- and dual-slopes) is possible by means of generalized pathlengths.

[0127] The methods and systems described herein are general and applicable to any geometry of diffusive medium and any distribution of the optical properties within the medium. It also has the advantage of not mixing data types with different features of spatial sensitivity. The sensitivity is defined for a single focal absorption change as the ratio of the equivalent absorption change to the true focal change. Equivalent absorption changes to focal absorption perturbations calculated with single-slope, feature regions of negative sensitivities, particularly close to the short distance detector when the single-slopes are calculated by using one source and two detectors.

[0128] By using a symmetric source-detector arrangement, the dual slope method also achieves a partial correction and cancellation of the negative sensitivity. This arrangement is similar to that used for a “self-calibrating” approach to absolute estimation of the optical properties by which it is possible to measure absolute optical properties of a diffusive medium without previous calibration. The self-calibrating method features high insensitivity to light source power and detectors sensitivity fluctuations as well high insensitivity to changes in the coupling between input source and the probed medium and medium and the detectors (usually trough optical fibers). The self-calibrating approach is also used for measuring dynamic changes of oxy-, deoxyhemoglobin, and oxygen saturation. In these measurements, the changes in DC slope originated from homogeneous changes in the absorption coefficient.

[0129] The methods and systems described herein result from a rigorous and general framework for this approach (not restricted to DC data only), that in principle could be used also for a more rigorous inversion problem (see for example equation (13) that correlates the change in slopes and the focal changes in the absorption coefficient).

[0130] For single-distance data, the sensitivity of DC is higher than that of the phase at shallow depths, while the opposite is true deeper in the medium. This is the case for single-distance, single-slope, and dual-slope data. However, DC data is characterized by a much greater SNR than phase data. Therefore, even though phase dual-slope has some of the most desirable features for imaging purposes, due to its low SNR it may not be always applicable to experimental data, depending on the nature (absorption contrast, size, and duration) of the target optical contrast. In this case, one can use the DC dual-slope method, which features a higher SNR and a better ratio of sensitivities between deep and shallow regions than DC single-distance data.

[0131] The dual-slope method requires a minimum of two sources and two detectors. There is no need for a short source-detector separation to suppress extracerebral hemodynamic contributions to signals collected at large source-detector distances. This is because the dual-slope method (especially with phase data) is intrinsically weakly sensitive to localized changes close to the surface. However, the dual-slope method can also be used with multiple source-detectors (see FIG. 1D) and therefore also benefit from the information of several single-distance data.

[0132] The theoretical results suggest that dual phase slope, despite the low signal-to-noise ratio (SNR), can be used for typical scenarios of near-infrared spectroscopy. A particular application is that of retrieval of the ratio of oxy- and deoxyhemoglobin oscillations in a typical protocol used in coherent hemodynamics spectroscopy as a way to obtain information on cerebral microvascular integrity and autoregulation. The method is also applicable to retrieving the ratio of oxy- and deoxyhemoglobin oscillations in a typical case of periodic brain activation. In both cases the dual-slope methods (particularly with phase data) yielded ratio of hemoglobin oscillations that were close to actual values in the brain. The dual-slope method could be used in other applications of near-infrared spectroscopy as a simpler alternative to more complex methods typical of diffuse optical tomography. As one example, the dual-slope method can be successfully applied to experimental data from the human head during a typical protocol of coherent hemodynamic spectroscopy.

[0133] Finally, a multi-distance formula used in the semi-infinite medium geometry uses both AC and phase slope to retrieve the absolute value of the absorption and the reduced scattering coefficients. The formula for the absorption coefficient can also be used for estimating equivalent absorption changes due to arbitrary focal changes in the medium. However, both phase slope and AC slope changes will contribute to the absorption changes and these two data types have very different spatial sensitivity, causing a poorer sensitivity to deeper regions than that obtained with dual phase slope data alone.

5. Conclusions

[0134] The sensitivity of single-distance data types and derived slopes to focal absorption changes demonstrates the advantages of a dual-slope method. While the spatial resolution in diffuse optics is intrinsically limited by the high ratio (a few hundreds) between source-detector separation and scattering mean free path, it is nevertheless possible to have a family of source-detector configurations that achieve a deeper and more localized region-of-sensitivity than that which is available using typical source-detector configurations of the type commonly used in diffuse optics.

[0135] In particular, a 55-millimeter linear array of two sources and two detectors promotes the ability to determine dual-slope intensity and phase data with maximal sensitivity at depths of five millimeters and eleven millimeters respectively under typical conditions of near-infrared spectroscopy of blood perfused tissue. This is a marked improvement over the sensitivity at depths of less than ten millimeters and less than five millimeters for single-distance intensity and phase data, respectively. This result points to more effective non-invasive measurements of brain and skeletal muscle under superficial layers of skin, skull, adipose tissue, etc.

[0136] In addition to a deeper sensitivity, the proposed dual-slope method features significant practical advantages that were previously demonstrated in a self-calibrating approach for absolute measurements of optical properties. These advantages include an insensitivity to instrumental drifts (source emission properties or detector sensitivity response), and variable opto-mechanical coupling between optical probe and tissue (as may result from subject movement or probe adjustment over time) that occur over a longer time scale than the time resolution of data collection. This is a crucially important feature to enhance the reliability of the data collected on living tissue over relatively long periods of time on moving subjects.

[0137] Having described the invention and a preferred embodiment thereof, what is claimed as new and secured by letters patent is: