1. Introduction

The first traceable attempt to use light to image biological tissues was made as early as the 1920’s [1], but at that time the problem of high scattering was not resolved. A new approach to optical image reconstruction as a form of inverse problem [2] emerged about 20 years ago [3,4]. Since then, a significant effort has been made to develop algorithms for image reconstruction (see a comprehensive review by S. Arridge [2]), and also multichannel imaging instruments (see a review by P. Rolfe [5] and references therein). However, practical achievements in three-dimensional tomography of adult human brain [6–11], i.e. the spatial localization of optical properties and their changes, have been significantly more modest than the results of purely spectroscopic [8, 12–17] and two-dimensional mapping [18–21] studies. The latter have exploited the ability of near infrared (NIR) spectroscopy to provide direct measurement of oxy- and deoxy-hemoglobin concentrations, as well as detection of the fast neuronal responses with millisecond time resolution [17]. The origin of the image reconstruction problem is two-fold. The first difficulty remains the same as in the 1920’s, i.e., the high degree of light scattering by tissues which, in spite of application of sophisticated algorithms and imaging devices, results in severe underdetermination of the inverse problem. The second complication is the general non-linearity of the inverse problem with respect to the optical properties of the medium [2].

In the case of measuring functional activity in the brain it is generally assumed that the non-linearity can be neglected since the expected changes in the optical properties are small [8]. This linearized approach is also known as the perturbation approach [2]. However, even after linearization, the underdetermination of the problem, i.e., the fact that the number of measurements (typically a few hundred) is significantly lower than the number of voxels (several thousand for the adult head), remains a significant problem. This problem is complicated by the fact that near-infrared sensors are much more sensitive to changes close to the locations of sources and detectors than in deep tissues such as the brain. Recently, approaches to this problem have been proposed in the form of a combination of NIR imaging with other structural imaging modalities to determine the geometric structure of deep-lying tissues, thus allowing restriction of the area for image reconstruction to the cortical region [9, 15, 22].

Any image reconstruction algorithm must include a solution to the forward problem of light transport [2]. Although the diffusion approximation [2,23] has been shown to be adequate for most tissues, thus providing a very efficient approach to solving the forward problem, a well-known obstacle to applying the diffusion or higher order deterministic approximations [2] in the adult brain imaging is the presence of low-scattering cerebro-spinal fluid (CSF) [24]. For this reason approaches based on a direct Monte Carlo (MC) simulation of photon migration have been developed by a number of groups [4,25]. Although the MC method is intrinsically slower than the numerical solution of partial differential equations (PDE), its robustness together with ongoing advances in computer technology offer a valuable alternative to PDE methods.

In a recent paper [9] the Massachusetts General Hospital Photon Migration group proposed a number of ways to improve the performance of diffuse optical tomography (DOT), namely the selection of particular wavelengths, implementation of overlapping detector channels, and the utilization of prior spatial information from magnetic resonance imaging (MRI). In this paper, we expand our previous work [15,16,20] and that of others [9,11,13,14] on simultaneous functional magnetic resonance imaging (fMRI) and NIR spectroscopy by developing a methodology that integrates simultaneously-acquired fMRI with DOT, incorporating several features which have also been outlined in simulation studies [11]. The key advances of our work are the design of an MRI-compatible optical probe with overlapping source/detector channels, Monte Carlo simulations of photon migration paths based on structural MR data, light sources at optimal wavelengths of 690 and 830 nm, and the restriction of the volume used for iterative reconstruction to brain tissue alone, based on structural MRI information. In this context the simultaneous measurement of optical and MRI data allows highly accurate and reproducible image co-registration, as well as providing a high-resolution subject-specific structural model for Monte Carlo simulations of photon migration paths.

2. Design of an MRI-compatible optical probe

The strong static magnetic field, rapidly switched gradients and radiofrequency (RF) pulses used in MRI measurements mean that the optical probe must be non-magnetic and, ideally, completely metal-free. The optical probe was designed with 16 pairs of 400-µm-diameter core plastic-clad multimode silica source fibers and 4 detector fiber bundles. The physical size of the 3T head-only MR scanner (Siemens Allegra) mandated the use of prisms with dielectric reflective layers for the fiber bundles (metal reflecting surfaces were found to produce significant artifacts in MR images). The optical fibers were “ferruled” using plastic tubing, and the frame of the optical probe was made polyurethane which has appropriate mechanical properties and produces insignificant artifacts in the MR images. The topology of the probe was designed so that the optical channels overlapped and the distribution of source-detector distances covers the optimal range – between approximately 20 and 30 mm, as shown in Fig. 1. An MRI visible marker was attached adjacent to each of the 16 optical source fibers and four detector fiber bundles so that accurate source and detector positions could be estimated from MR images. The thickness of the probe was less than 20 mm to ensure that it could be placed comfortably between the back of the head and the bottom of the MRI birdcage head coil. For the experiment, the center of the optical probe was placed over the primary visual cortex of the subject. A pair of non-magnetic goggles (Resonance Technology, Inc.) with LCD screens were placed in front of the subject’s eyes inside the birdcage head coil. The visual stimulation paradigm consisted of five blocks each with 28.8 s fixation followed by 19.2 s of a black-and-white checkerboard pattern flashing on for 50 ms and off for 1.95 s.

 

Fig. 1. (a) Schematic of the fMRI-compatible optical probe. (b) Optode positions reconstructed from MR images, superimposed on a surface rendering of a 3D MRI data set.

Download Full Size | PDF

3. Data acquisition

The optical signals were recorded using a near infrared spectrometer (Imagent, ISS, Champaign, IL). Data acquisition was synchronized with the fMRI measurements using the TTL-trigger signal from the MR scanner, which also triggered the beginning of the visual stimulation paradigm. The optical sources were laser diodes (690 and 830 nm) which were amplitude modulated at 150 MHz and time-multiplexed. Light reaching the detectors was amplified by photomultiplier tubes and consequently converted into AC, DC, and phase signals for each of the source-detector combinations, or channels, at each wavelength.

The fMRI protocol started with a rapid T₁-weighted spin-echo (240×240 mm field-of-view, 512×512 data matrix, 16 slices, 4 mm slice thickness, 0.4 mm inter-slice gap, aligned along the anterior/posterior commisure). These images were used for the initial co-registration with the functional scans, which were acquired using an echo planar imaging (EPI) pulse sequence (slice orientation and dimensions as above, 64×64 data matrix, tip angle 60°, TE 25 ms, TR 2000 ms). After the functional data have been collected, high-resolution (0.5×0.5×1.0 mm³) T₁-weighted spin-echo images, termed “probe localizers”, were acquired over a small volume containing the optical probe to enable automatic recovery of the positions of the optodes. Finally, a magnetization-prepared rapid gradient-echo (MPRAGE) 3-D image (spatial resolution 0.94×0.94×1.2 mm³) of the full-head was acquired for co-registration of the optical and fMRI signals.

4. Data processing and optical image reconstruction

Automatic estimation of source- and detector-positions from the MR images, incorporating deformation of the optical probe, was performed using custom-developed image processing algorithms written in MATLAB (Mathworks, Natick, NJ). The process can be summarized as computing the rotation and translation of the reference coordinate systems of the MPRAGE image with respect to the probe-localizer, determination of the angular position of the cylindrical markers by computing the rotation of the coordinate system due to rotation of the optical probe, and finally correcting for any misalignment from deformation of the optical probe. The accuracy of the optode localization was comparable with the resolution of the structural MR images (~1 mm). An example of the estimated optode positions is shown in Fig. 1(b).

The MPRAGE images were also used, after segmentation, as the template for the Monte Carlo simulations of light transport through the subject’s head [25]. Although Monte Carlo simulations are computationally expensive, they avoid known problems related to non-scattering regions, such as cerebrospinal fluid, and to the breakdown of the diffusion approximation near boundaries. According to the Monte Carlo model the complex fluence from source s, modulated at frequency ω and measured by detector d, can be expressed in terms of the sum of the complex amplitudes of light traveling along all possible paths p in the scattering medium:

where, in a given voxel n, $l_n^p$ is the total pathlength of the light, which can be larger than the voxel size since the path may pass through the voxel many times, and k_n is the complex wavevector of the medium given by, in terms of ${{\mu }_{\mathrm{a}}}^{\mathrm{n}}$ the absorption coefficient, k_n =-${\mathit{\mu }}_a^n$ +iω/c. Assuming that cerebral activity causes the local absorption coefficient to change from its baseline value, $\overline{\mu }$ⁿ_a, by a small value ${\delta \mathit{\mu }}_a^n$ , the change in the fluence can be easily obtained by applying a linear approximation to Eq. (1):

where k̄_n is the baseline wavevector, index m is related to summations over only the activated voxels, and index p is related to summations along those paths which pass from the source to the detector through the activated voxels. In order to compute the fluence perturbation using Eq. (2) one needs to simulate baseline conditions only for those paths from the source (and not from the detector as in Ref. [26]) to the particular voxel, and to determine only those paths which pass through the imaging volume to the detector, a process which reduces significantly the required computational memory and time. One can also show that Eq. (2) is consistent with the Born approximation under conditions satisfying the diffusion equation. It should be noted that the expression in square brackets in Eq. (2) represents the complex light bundle, K_m , or the sensitivity point-spread function of a source-detector pair [2].

A commonly-used and well-documented Monte Carlo software package [25] was modified to store the trajectory of each detected “light pathway”, and was improved in terms of the required computational resources for frequency-domain measurements. In order to validate this approach, we have numerically verified that, in the case of a homogeneous semi-infinite highly-scattering medium, the complex light bundle given by K_m corresponds to that computed using an analytical solution to the diffusion equation [27]. We found that 10⁷ simulated light paths provides a signal-to-noise (S/N) of ~20:1 near the pulse peak in the human head, representing a good compromise between S/N and required computational resources. The pulses were computed within 10 ns time with 50 ps temporal resolution. The noise was approximately stationary between 0.2 and 1 ns. The S/N was estimated assuming an averaging time of 0.2 ns near the pulse peak and a bandwidth of 10 GHz. The simulations were performed over the entire volume of the head and the required computation time for each source-detector pair was 400 minutes on a 2.5 GHz Opteron (64-bit) CPU.

In order to compare the experimental measurements δV with predicted values of δU _sd in Eq. (2), δV must be normalized by the baseline value V̄, and δU _sd by the baseline simulated value Ū_sd computed using Eq. (1). This results in

Complex data acquired for each source-detector combination as a function of time are termed a measurement. For M measurements over an imaging volume of N voxels, Eq. (3) can be written in a generalized matrix form:

The M-by-N coefficient matrix A consists of the bundle K_m evaluated for each of N voxels and each of M measurements. The M-by-1 measurement vector B consists of the detected optical signal. The solution X is a 1-by-N vector that contains the changes in the absorption coefficient of all the voxels within the imaging volume. Eq. (4) is typically under-determined and ill-posed. For the spatial resolution of the optical image reconstruction used in this study (4×4×4 mm³), the number of measurements (M=64) is much less than the number of voxels (N=10³~10⁴). Applying physiological and spatial a priori constraints can effectively reduce the number of unknowns [9]. Here, we have constrained the solution to the “optical sensitivity region” in the brain. In order to determine such a sensitivity region we computed signal changes due to perturbation of the absorption coefficient in each voxel in the brain and compared the result with the experimentally assessed measurement error caused by the noise. Assuming an “activation level” of 5% of the baseline absorption (a conservative estimate based on previous NIRS studies of functional activation in humans, such as in Ref. [16]), we considered a voxel to be detectable by a particular measurement channel if it produced a signal change larger than the corresponding standard error. This error was computed as the standard deviation of the raw data noise (i.e., the optical signal changes at baseline conditions; statistical assessment using the Matlab Statistical Toolbox has revealed that this noise is Gaussian) divided by the square root of the number of blocks in the brain activation paradigm, in this case five. Typically, this reduces the number of voxels to be reconstructed to ~600. A number of image reconstruction algorithms were assessed, including matrix left division, singular value decomposition (SVD), algebraic reconstruction techniques (ART), and simultaneous iterative reconstruction technique (SIRT) for both simulated and experimental data. Among the above algorithms tested, SIRT gave the best result in terms of the quality of the reconstructed optical images.

5. Simulation results

Simulations and image reconstructions were conducted based on photon migration in both homogeneous and heterogeneous models of the head, the latter generated by segmenting a 3-D anatomical MR image of the head into scalp/skull, cerebrospinal fluid (CSF), gray matter and white matter, and assigning baseline optical properties to each tissue type based on values in the literature [2,25]: the values of µ_a (mm^-1) were 0.014 for scalp and skull, 0.004 for CSF, and 0.01 for both gray and white matter, corresponding values of µ_s (mm^-1) were 0.86, 0.1 and 1.11, the refractive index was set to 1.4 for all tissues, and since no accurate literature values of the scattering anisotropy factor g were available, it was assumed that g=0, corresponding to an almost homogeneous distribution of the scattering angle.

In order to simulate hemodynamic activation within the head, a sphere with a diameter of 20 mm was embedded into a digital head model derived from MR data. The optical properties of the sphere were the same as gray matter, with a 10% increase in absorption coefficient when “activated”. The sphere was positioned within the brain close to the surface of the visual cortex in the left hemisphere, approximately 5 mm to the left of the center of the optical probe (see Fig. 2(a)). The results of the image reconstruction from simulated measurements using SIRT are shown in Figs. 2(b), 2(c), and 2(d). In Fig. 2 one can see that in the superficial cortical region reconstruction is reasonably accurate. Deeper in the brain the accuracy decreases significantly as a result of the lower light density. The SIRT algorithm is well-known to assign the largest values where the sensitivity is largest, i.e., where the bundle is most dense, which is the surface of the brain.

 

Fig. 2. a) Simulated activation (modeled as a change in absorption coefficient) compared with the (b) reconstructed result. One-dimensional projections of the simulated and reconstructed absorption perturbations are plotted as functions of distance from the center of the cross-hairs along the X (c) and Y (d) coordinates.

Download Full Size | PDF

6. Human experimental results

The identical optical image reconstruction procedure was performed for data acquired on human subjects with the visual stimulation paradigm described previously. Using Lambert-Beer’s law, the changes in absorption coefficient were converted into changes in oxy- and deoxyhemoglobin concentrations. For each voxel, the correlation coefficient between changes in oxy- and deoxyhemoglobin concentration and the stimulation paradigm, delayed by six seconds, was calculated. The spatial hemodynamic responses of the oxy- and deoxyhemoglobin are compared with the fMRI results in Figs. 3(a)–3(c). The fMRI data were processed using the Automated functional neuroimaging (AFNI, Medical College of Wisconsin) and statistical parametric modeling (SPM) programmes as follows: the data were motion corrected (AFNI), detrended (AFNI), and then the correlation and percentage change maps obtained (AFNI). Voxels with BOLD signal change higher than 8% were removed as they represent activation in large veins. The BOLD image in Fig. 3 shows the percentage change thresholded by a correlation value of 0.5, which corresponds to a p value less than 0.01 [28]. The functional maps were registered to the 3D MPRAGE dataset (SPM) via intermediate registration with the T1-weighted images.

 

Fig. 3. (a–c) Spatial maps of the correlation coefficient, shown as the colour bar on the right. (a) Spatial hemodynamic response derived from the BOLD MRI signal, (b) spatial hemodynamic response corresponding to the change in deoxyhemoglobin concentration derived from the reconstructed optical data, (c) spatial hemodynamic response corresponding to the change in oxyhemoglobin concentration derived from the reconstructed optical data. (d–e) Temporal changes in oxy- and deoxyhemoglobin concentrations derived from the optical data. (d) Temporal hemodynamic response of corresponding to an average of all voxels with a correlation coefficient greater than 0.5, (e) response from voxels within the sensitivity region but with a correlation coefficient less than 0.5, and (f) response from the single voxel showing the maximum changes in the oxy- and deoxyhemoglobin concentrations.

Download Full Size | PDF

Reasonable consistency can be seen between the spatial hemodynamic response of the optical signals and the fMRI signal. Temporal changes in oxy- and deoxyhemoglobin concentrations were averaged over “activated” voxels that had a correlation coefficient greater than 0.5. The changes in oxy- and deoxyhemoglobin concentration demonstrated very similar temporal patterns, as shown in Fig. 3(d). Data for non-activated regions are also shown in Fig. 3(e).

To conclude, a complete functional imaging technique of simultaneous MRI and DOT has been implemented and tested in vivo. To the best of our knowledge, this is the first in vivo optical imaging result obtained using a Monte Carlo approach of a patient-specific head model, iterative reconstruction, and direct comparison with simultaneously acquired fMRI data. Although DOT images are inferior to those acquired using MRI in terms of spatial resolution, the strength of DOT is in providing biochemically specific time-resolved and spatially localized information on cerebral processes. In addition to humans, this technique can be applied readily to small animal imaging, where the small size of the brain requires source-detector distances less than 5 mm, at which the diffusion approximation gives poor results and the Monte Carlo approach is necessary.