Three-dimensional optical tomography of hemodynamics in the human head

We report on the first three-dimensional, volumetric, tomographic localization of vascular reactivity in the brain. To this end we developed a model-based iterative image reconstruction scheme that employs adjoint differentiation methods to minimize the difference between measured and predicted data. The necessary human-head geometry and optode locations were determined with a photogrammetric method. To illustrate the performance of the technique, the three-dimensional distribution of changes in the concentration of oxyhemoglobin, deoxyhemoglobin, and total hemoglobin during a Valsalva maneuver were visualized. The observed results are consistent with previously reported effects concerning optical responses to hemodynamic perturbations. ©2001 Optical Society of America OCIS codes: (170.6960) Tomography; (170.3880) Medical and biological imaging; (170.0110) Imaging systems; (170.1470) Blood/tissue constituent monitoring ______________________________________________________________________________________________ References and links 1. D. A. Benaron, S. R. Hintz, A. Villringer, D. Boas, A. Kleinschmidt, J. Frahm, C. Hirth, H. Obrig, J. C. Van Houten, E. L. Kermit, W. Cheong, D. K. Stevenson, “Noninvasive functional imaging of human brain using light,” J. Cerebral Blood Flow and Metabolism 20, 469-477 (2000). 2. Y. Hoshi, I. Oda, Y. Wada, Y. Ito, Y. Yamashita, M. Oda, K. Ohta, Y. Yamada, M. Tamura, “Visuospatial imagery is a fruitful strategy for the digit span backward task: a study with near-infrared optical tomography,” Cognitive Brain Research 9, 339-342 (2000). 3. E. Watanabe, A. Maki, F. Kawaguchi, Y. Yamashita, H. Koizumi, Y. Mayanagi, “Noninvasive cerebral blood volume measurement during seizures using multichannel near infrared spectroscopic topography,” J. Biomed. Opt. 5, 287-290 (2000). 4. M. Franceschini, V. Toronov, M. E. Filiaci, E. Gratton, S. Fantini, “On-line optical imaging of the human brain with 160-ms temporal resolution,” Optics Express 6, 49-57 (2000), http://www.opticsexpress.org/oearchive/source/18957.htm 5. S. Fantini, D. Huebert, M. A. Franceschini, E. Gratton, W. Rosenfeld, P. G. Stubblefield, D. Maulik, and M. Stankovic, “Non-invasive optical monitoring of the newborn piglet brain using continuous-wave and frequency-domain spectroscopy,” Phys. Med. Biol. 44, 1543-1563 (1999). 6. M. R. Stankovic , D. Maulik, W. Rosenfeld, P. G. Stubblefield, A. D. Kofinas, S. Drexler, R. Nair, M. A. Franceschini, D. Hueber, E. Gratton, and S. Fantini,, “Real-time optical imaging of experimental brain ischemia and hemorrhage in neonatal piglets,” J. Perinat. Med. 27, 279-286 (1999). 7. H. Koizumi, Y. Yamashita, A. Maki, T. Yamamoto, Y. Ito, H. Itagaki, and R. Kennan, “Higher-Order Brain Function Analysis by trans-cranial dynamic near-infrared spectroscopy imaging,” Journal of Biomedical Opt. 4, 403-413 (1999). 8. A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. 20, 435-442 (1997). 9. M. Tamura, Y. Hoshi, and F. Okada, "Localized near-infrared spectroscopy and functional optical imaging of brain activity," Philosophical Transactions of the Royal Society of London Series B: Biological Sciences. 352(1354), 737-42 (1997). (C) 2001 OSA 10 September 2001 / Vol. 9, No. 6 / OPTICS EXPRESS 272 #34858 $15.00 US Received July 27, 2001; Revised August 29, 2001 10. A. Kleinschmidt, H. Obrig, M. Requardt, K. Merboldt, U. Dirnagl, A. Villringer, and J. Frahm, “Simultaneous recording of cerebral blood oxygenation changes during human brain activation by magnetic resonance imaging and near infrared spectroscopy,” J. of Cerebral blood flow and metabolism 16, 817-826 (1996). 11. H. Liu, B. Chance, A. H. Hielscher, S. L. Jacques, F. K. Tittel, “Influence of blood vessels on the measurement of hemoglobin oxygenation as determined by time-resolved reflectance spectroscopy,” Med. Phys. 22, 1209-1217 (1995). 12. G. Gratton and M. Fabiani, “Dynamic brain imaging: Event-related optical signal (EROS) measures of the time course and localization of cognitive-related activity,” Psychonomic Bulletin and Review 5, 535-563 (1995). 13. A. Maki, Y. Yamashita, Y. Ito, E. Watanabe, Y. Mayanagi, H. Koizumi, “Spatial and temporal analysis of human motor activity using non-invasive NIR topography,” J. Med. Phys. 22, 1997-2005 (1995). 14. Y. Hoshi and M. Tamura, “Detection of dynamic changes in cerebral oxygenation coupled to neuronal function during mental work in man,” Neuroscience Letters, 150, 5-8 (1993). 15. S. R. Hintz, W. F. Cheong, J. P. van Houten, D. K. Stevenson, D. A. Benaron, “Bedside imaging of intracranial hemorrhage in the neonate using light: comparison with ultrasound, computed tomography, and magnetic resonance imaging,“ Pediatr. Res. 45, 54-59 (1999). 16. A. H. Hielscher, A. D. Klose, K. M. Hanson, “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imag. 18, 262-271 (1999). 17. A. D. Klose, and A. H. Hielscher, “Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer,” Medical Physics 26, 1698-1707 (1999). 18. R. Roy and E. M. Sevick-Muraca, “Truncated Newton's optimization scheme for absorption and fluorescence optical tomography: Part I Theory and formulation,” Opt. Express 4, 353-371 (1999). http://www.opticsexpress.org/oearchive/source/9268.htm 19. S. R. Arridge, M. Schweiger, “A gradient-based optimization scheme for optical tomography,” Opt. Express 2, 213-226 (1998). http://www.opticsexpress.org/oearchive/source/10447.htm 20. M. Schweiger, S. R. Arridge, M. Hiraoka, D. T. Delpy, “The finite element model for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22, 1779-92 (1995). 21. S. R. Arridge, M. Schweiger, M. Hiraoka, D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20, 299-309 (1993). 22. H. Jiang, K. D. Paulsen, and U. L. Osterberg, “Optical image reconstruction using DC data: simulations and experiments,” Phys. Med. Biol. 41, 1483-1498 (1996). 23. G. Abdoulaev, A. Varone and G. Zanetti, “An object oriented flow solver for the CRS4 virtual vascular project,” in Proc. Of the 3 European Conf. on Numerical Mathematics and Advanced Applications, P. Neittaanmaki, T. Tiihonen and P. Tarvainen (eds.) (World Scientific, Singapore), 407-415 (2000). 24. J. H. Bramble, J. E. Pasciak, and J. Xu, “Parallel multilevel preconditioners,” Math. Comp. 55, 1-22 (1990). 25. Y. Pei, H. L. Graber, R. L. Barbour, “Influence of systematic errors in reference states on image quality and on stability of derived information for DC optical imaging,” Appl. Opt. 40 (2001) in press. 26. D. A. Boas, T. Gaudette, G. Strangman, X. Cheng, J. J. Marota, J. B. Mandeville, “The accuracy of near infrared spectroscopy and imaging during focal changes in cerebral hemodynamics,” Neuroimage 13, 7690 (2001). 27. A. J. Davies, D. B. Christianson, L. C. W. Dixon, R. Roy, P. van der Zee, “Reverse differentiation and the inverse diffusion problem,” Advances in Engineering Software 28, 217-221 (1997). 28. D.B. Christianson , A. J. Davies, L. C. W. Dixon, R. Roy, P. Van der zee, “Giving Reverse Differentiation a Helping Hand,” Optimization Methods and Software 8, 53-67 (1997). 29. H. M. Karara, Handbook of Non-topographic Photogrammetry (American Society of Photogrammetry, Falls Church, VA, 1989). 30. E. F. Church, Elements of Photogrammetry (Syracuse University Press, New York, 1948). 31. G. Deng, W. Falg, “An evaluation of an off-the-shelf digital close-range photogrammetric software package,” Photogrammetric Engineering and Remote Sensing 67, 227-233 (2001). 32. R. H. Bartels, J. C. Beatty, and B. A. Barsky, An Introduction to Splines for Use in Computer Graphics and Geometric Modeling (Morgan Kaufman Publishers, 1987). 33. C. H. Schmitz, M. Löcker, J. Lasker, A. H. Hielscher, R. L. Barbour, "Performance characteristics of a silicon-photodiode-based instrument for fast functional optical tomography," in Optical Tomography and Spectroscopy of Tissue IV, B. Chance, R. R. Alfano, B. Tromberg, M. Tamura, E. M. Sevick-Muraca, Proc. SPIE 4250, 171-179 (2001). 34. C. H. Schmitz, H. L. Graber, H. Luo, I. Arif, J. Hira, Y. Pei, A.Y. Bluestone, S. Zhong, R. Andronica, I. Soller, N. Ramirex, S. L. S. Barbour, R. L. Barbour, “Instrumentation and calibration protocol for imaging dynamic features in dense-scattering media by optical tomography,” Appl. Opt. 39, 6466-6486 (2000). 35. F. P. Tiecks, C. Douville, S. Byrd, A. M. Lam, D. W. Newell, “Evaluation of impaired cerebral autoregulation by the Valsalva Maneuver,” Stroke 27, 1177-1182 (1996). 36. X. Cheng, D.A.Boas, “Systematic diffuse optical image errors resulting from uncertainty in the background optical properties,” Optics Express 4, 299-307 (1999), http://www.opticsexpress.org/oearchive/source/9108.htm (C) 2001 OSA 10 September 2001 / Vol. 9, No. 6 / OPTICS EXPRESS 273 #34858 $15.00 US Received July 27, 2001; Revised August 29, 2001 37. S. Wray, M. Cope, D. T. Delpy, “Characteristics of the near infrared absorption spectra of cytochrome aa3 and hemoglobin for the noninvasive monitoring of cerebral oxygenation,“ Biochim Biophys Acta 933, 184192 (1988). 38. T. O. McBride, B. W. Pogue, U. L. Osterberg, K. D. Paulsen, “Separation of absorption and scattering heterogeneities in NIR tomographic imaging of tissues," in Biomedical Topical Meetings, OSA Technical digest (Optical Society of America, Washington DC, 2000), pp. 339-341. 39. Y. Pei, H. L. Graber, R. L. Barbour, “Normalized-constraint algorithm for minimizing inter-parameter crosstalk in DC optical tomography,” Opt. Express 9, 97-109 (2001), http://www.opticsexpress.org/oearchive/source/33900.htm 40. Y. Hoshi, N. Kobayashi, M. Tamura, “Interpretation of near-infrared spectroscopy signals: a study with a newly developed perfused rat brain model,” J. Appl. Physiol. 90, 1657-1662 (2001). 41. F. P. Tiecks, A. M.

[1] studied physiological changes in brain oxygenation in male adults during mixed motor and sensory cortex activation.Optical imaging was performed during a sequential thumb-tofinger apposition task.The baseline brain oxygenation, being stable over time, was subtracted from an activated state (motor task) to unmask a signal.Y. Hoshi et al. [2] obtained quantitative images of hemoglobin concentration changes associated with neuronal activation in the human brain during a forward (DF) and backward (DB) digit span task, which assesses verbal working memory.Topographic subtraction images of changes in oxyhemoglobin were generated between a task and its previous resting state.These images were then superimposed on 3-D structural MRI maps of the head as a means of ascertaining spatial localization.E. Watanabe et al. [3] have developed an optical system for obtaining 2dimensional maps localizing epileptic foci.This method determines locations of blood flow increase during epileptic episodes.Franceschini et al. [4] have reported on imaging arterial pulsation and motor activation in healthy human subjects.These and other works have established that near-infrared (NIR) light can be used to probe the brain for changes in blood oxygenation and blood volume.
In most of the currently employed protocols, either a single source and a single detector is used to perform point measurements [10,14], or an array of sources and detectors with a fixed geometry is employed to obtain topographic maps [1][2][3][4]6,7,13,15].These maps are typically generated by back-projecting the changes in the absorption coefficients along the straight-line paths between sources and detectors [2,4].The resulting surface maps project the cortical response during various protocols together with the superficial vascular changes onto a twodimensional surface map.
The goal of the presented work has been to go beyond topographic maps and obtain threedimensional, volumetric reconstructions of the changes in optical properties inside a human head.To this end we have developed a model-based iterative image reconstruction (MOBIIR) scheme that uses a finite-element formulation of the three-dimensional, time-independent diffusion equation as a forward model.To apply MOBIIR schemes to brain imaging, the surface geometry and the locations of sources and detectors on the surface of the head have to be determined.Here we describe our approach to image recovery, as well as a convenient way to determine the surface geometry of irregular structures.To illustrate the performance of this approach we conducted measurements on the human forehead during a Valsalva maneuver, and generated volumetric reconstructions that display the hemodynamics during the maneuver.

Image Reconstruction Algorithm
In this work the 3-dimensional reconstruction of the optical properties in the human head was achieved using a MOBIIR scheme [16][17][18][19].In general, a MOBIIR scheme comprises three major parts: (1) a forward model that predicts the detector readings based on an assumed spatial distribution of optical properties, (2) an objective function Φ that compares predicted with measured signals, and (3) a scheme that updates the assumed optical parameters for subsequent forward calculations.
The forward model that we employ is the time-independent diffusion equation with Robin boundary conditions along all exterior surfaces.The finite element forward model is similar to works that have already been described by other groups [18,[20][21][22].The domain Ω is divided into P elements, joined at N vertex nodes.The solution is approximated by a piecewise linear function spanned by a set of linear basis functions .
ϕ The solution is obtained by sparse matrix inversion using the preconditioned conjugate gradient method [23].Two added features allow for mesh refinement and multigrid capabilities [24].The discretized forward model can be written as c A S φ = where A is the finite element stiffness matrix of the diffusion equation forward model.Each element (el) of the tetrahedral mesh, with nodes i, j has a corresponding 4 4 × A e element stiffness matrix with entries defined by: where the diffusion coefficient 1/ 3( ') , µ a is the absorption coefficient, and ' s µ is the reduced scattering coefficient.Our formulation of the MOBIIR scheme requires that predicted detector readings be compared to actual measurements.To this end it is necessary to define an objective function that determines the goodness-of-fit between measured data, M, and predicted detector data, P.Many groups have consistently used the least-square objective function 2 , , 2 , ( ) where s, d refer to sources and detectors, respectively and σ is the standard deviation [16][17][18][19].
However, to use Eq. 2 absolute measurements of M are necessary, which are often difficult to obtain.More commonly, researchers in near-infrared imaging generate data via a differencemeasurement approach.In a difference-measurement approach one is unable to determine the absolute detector readings for a single time point and hence must compare the change in detector readings between two states.This approach has two main advantages: first, it is less sensitive to boundary effects, and second, it is less sensitive to the initial guess chosen for the background medium [25].Its major disadvantage is that one cannot determine the true distribution of optical properties, only the change in the µ a , µ s ', or D from a given baseline.
However, many groups have used similar difference-measurement approaches for localizing brain activity and for determining general trends in the oxygenation state [1,8,26].
To deal with difference data, we have adapted the approach suggested by Y. Pei et al. [25], and modified the objective function defined in Eq. (2) to yield where M pert is obtained during the perturbation, and M ref is a predetermined reference state.
Here the relative changes .This gradient defines a direction of steepest descent, which is then used in a multidimensional conjugate gradient descent algorithm that minimizes the objective function [16,17].
The calculation of the gradient was performed using the technique of adjoint differentiation [16, 17, 18, 27, and 28].In particular we adapted an approach described by Roy et al. [18].In our case this leads to the following formulation.Using the chain rule one can write the derivative of the objective function with respect to µ a for all nodes r as: where ( ) is the intensity at boundary node x d .The second term in parenthesis can be written as the dot product of two vectors where the first vector P r has a one at node r and zeros elsewhere: Combining Eq. ( 4) and (5) and rearranging the sum over sources and detectors one arrives at: The derivative of the last term in the outer parenthesis can be determined by taking the derivative of forward model A S ϕ = : Rewriting Eq. ( 7): Inserting Eq. ( 8) in Eq. ( 6) one arrives at: Defining the adjoint variables with * and noting that the sum is non-zero only at the detectors: After multiplying both sides by A, Eq. ( 10) is solved in a similar fashion to the forwardmodel equation.Combining Eq. ( 10) and Eq. ( 9) the derivative becomes: Using the fact that * * * , , , ( )

∑ ∑ ∑
, Eq. ( 11) can be rewritten as: One can interchange the sum over the source since , Pulling out the common term , This expression is computationally more efficient than Eq. ( 12) since the multiplication with where i ϕ are the basis functions, x k are the quadrature nodes, and I are the linear interpolation functions for the values at the nodes given by 4 So Combining Eq. ( 17) and Eq. ( 14) one arrives at the gradient of the objective function with respect to all optical properties µ a , ( ) . A similar mathematical formulation is derived to calculate the derivative of the objective function with respect to D. The computational burden is on the order of a forward calculation since the first term in Eq. ( 14) is calculated via Eq.( 10) (i.e., a "forward-like" calculation) and the second term is derived analytically, which is computationally efficient.Having generated gradients with respect to both sets of optical properties, these individual gradients are combined into an overall gradient of length 2n, where n is the number of nodes in the mesh (i.e. the number of optical parameters of one type).This procedure allows for the simultaneous update of both µ a and D. The aforementioned objective function and gradient algorithm is incorporated in a multidimensional conjugate gradient descent algorithm where typically 15-30 iterations with different gradients are necessary before convergence is achieved.

Determination of surface coordinates
For model-based optical tomographic imaging reconstructions one has to determine the surface geometry and optode locations.To achieve this we employed the techniques of photogrammetry [29,30,31].To determine the 3-D coordinates of the surface of the forehead as well as all optode locations we placed circular reflective markers on the head, as seen in Figure 1.These markers, because of their circular shape, have the advantage of allowing sub-pixel precision when used in combination with the centroid-locating algorithm.At the completion of the experimental data acquisition the pressure of the optodes against the skin left a mark that was easily identifiable and allowed for the accurate positioning of the markers on the forehead.After these 15 markers (corresponding to each optode) were positioned, the rest of the forehead was covered with markers that were used to determine the surface geometry of the entire forehead.Because the surface of the forehead can be approximated as a smooth contour it is acceptable to place a relatively sparse distribution of points over the surface of the head and use non-uniform B-splines (Nurbs) [32] to connect them.Other approaches, such as the use of holographic imaging or laser surface scanning devices are also feasible.However, the technical complexity paired with their high costs make them less practical and appealing in a day-to-day clinical setting.
Once all markers were placed on the surface of the head, the area of interest was photographed from different angles with a digital camera (SONY Mavika FD-90).These images of the head were then transferred to a desktop computer for analysis with the commercially available Photogrammetry software package (Photomodeler, Eos Systems Inc., Vancouver, Canada).After the targets in the photographs were marked and referenced we processed the images and computed the 3D coordinates of all referenced points.Using standard geometric models, we determined that the resulting coordinates on average have an accuracy of 0.25 mm with a standard deviation of 0.05 mm.The accuracy is dependent on the number of points referenced and on the maximum angle between the multiple camera positions, and a higher accuracy can be achieved using more reference points and more photographs taken from different angles.The reflective marker technique is suitable for the forehead and other body parts with little or no hair.In other circumstance, such as the hair-covered areas of the head, we have developed other referencing techniques that allow for the determination of the surface coordinates [33].For example, one can form a thermoplastic model with imbedded optodes to the contour of the head and then photograph the resulting mold, which is free of hair.

Mesh Generation
Having generated the surface coordinates of the human forehead, the boundary of our domain, it was necessary to build the corresponding volume.To this end, the coordinates of the surface points obtained using the photogrammetric approach were exported to a volumetric mesh generator.For this procedure, we used the CAD style volume generator GID software package (CIMNE International Center for Numerical Methods, Barcelona, Spain; http://gid.cimne.upc.es/intro/index.html).This package allowed us to manually extend the surface in the z direction by 4.0 cm and build the corresponding volume mesh, which is required for the reconstruction.We limited the volume to a 4.0 cm zone beneath the surface of interest, and not to the entire head, because we found that larger volumes are not necessary to consider, since their influence on the reconstruction results is negligible.This volume includes all areas with a minimum fluence value seven orders of magnitude smaller than the source strength.

Instrumentation and data acquisition
Measurements on the human forehead are performed with the dynamic near-infrared optical tomography (DYNOT) instrument, recently developed by Schmitz et al. [33,34].This system employs frequency-encoded multi-wavelength DC illumination, a time-multiplexed source, and parallel multi-channel detection together with fast gain switching (dynamic range of 180 dB).Fully configured, the basic unit can provide four wavelengths of simultaneous illumination at each illuminating site and can collect data from 32 channels in parallel at a source switching rate of ~75 Hz.For the current study, light produced by laser diodes operating at wavelengths of 760 nm and 830 nm was employed.Source light was intensity modulated at 5 KHz and 7 KHz, respectively, and optically coupled into transmitting fibers so as to provide for simultaneous dual-wavelength illumination at each site.This configuration allows for the discrimination between variations in oxyhemoglobin and deoxyhemoglobin levels as described in section 2.5.For the measurements presented in this work a three-tiered band was used to secure 3 sets of 5 optodes to the left forehead (Fig. 1).Four of the 15 optodes (labeled 3, 6, 10, 13 in Fig. 1) consisted of a co-located source and detector, while the remaining optodes were used only as detector sites (labeled 1-15 in Fig. 1).Tomographic measurements, involving all 4 15 × = 60 source-detector pairs, were performed at a frame rate of ~3 Hz throughout the acquisition period lasting approximately ten minutes.

Experimental Protocol
The experiment was designed to look at functional hemodynamic changes in the forehead of a single human subject induced by a Valsalva maneuver.During the Valsalva maneuver a forced expiration against a closed glottis demonstrates the effects of changes in intrathoracic pressure on blood pressure, and the brain's autoregulatory response to decreased vascular perfusion pressure in cerebral vessels [35].For the measurement the subject was placed in the supine position.Three epochs consisting of Valsalva maneuvers with one-minute rest periods interspersed were performed.
Assuming that the primary influences on the changes in the absorption coefficients at each wavelength are a linear combination of oxyhemoglobin and deoxyhemoglobin [36,37], one arrives at where λ indicates the wavelength, and ε Hb , ε Hb02 are the known extinction coefficients for deoxyhemoglobin and oxyhemoglobin at the given wavelengths [37], respectively.By simultaneously solving a set of algebraic equations at the disparate wavelengths 1 2,..., ( , ) it is possible to calculate the true concentrations of the, n, chromophores of interest.For the case of two chromophores consisting of oxyhemoglobin and deoxyhemoglobin the set of equations is: The values of a µ ∆ in Eq. (19-20), at each of the two wavelengths, are the calculated changes in the absorption coefficients at each node of the mesh determined using the reconstruction algorithm.Since the reconstruction algorithm determines both ∆µ a and ∆D, it explicitly accounts for the path length between given source detector positions and reconstructs the predicted values of a µ ∆ at each node.By solving Eq. ( 19) and Eq. ( 20), the changes in oxyhemoglobin, , can be determined.It should be noted that studies by others [19,20,21,25,38] and us [39] have shown that cross talk between the µ a and D reconstruction can occur.Hence, some of the oxyhemoglobin and deoxyhemoglobin effects may be attributed to scatter changes.Future works will have to address this issue either by modifying the reconstruction scheme or using different types of measurement strategies.

RESULTS
Traces of the measured output produced by the DYNOT system (section 2.1) at λ = 760 nm for the three epochs is shown in Fig. 2.This figure displays 15 traces of detector readings during three consecutive Valsalva maneuvers given a source at position number 3. Each trace was normalized by calculating the mean of the measured intensity during the first 30 seconds (rest period), and then dividing all intensities in the trace by that mean.The flat red trace corresponds to the detector co-located with the source at position number 3. Furthermore, the original data was filtered to enhance the signal to noise ratio, by using the median filter 1 2  For our reconstructions, the first of the three epochs was chosen, which is displayed in Fig. 3.At 15 seconds the Valsalva maneuver was begun.Initially a small decrease in signal for most detectors was observed, which returned towards baseline within approximately 10 seconds.Between the time points indicated by T1 and T6 the signal for all detectors decreased with a constant slope.At time point T6 a marked drop of signal intensity occurred for most detectors.A minimum was reached at about 65 seconds (T10) at which point the Valsalva maneuver was released.Subsequently, a rapid return towards base line of all traces occurred for 5 seconds, followed by a slower tapering towards the initial rest-period state.
For the volumetric difference reconstructions we used the ratio of the measured intensity during a time point (black arrows labeled T 1 -T 10 ) with respect to the measured intensity during the rest-state, red arrow (Fig. 3).These black arrows represent a snapshot (time point) of measured intensity for all source/detector combinations during consecutive instances of the Valsalva maneuver.For each snapshot a volumetric reconstruction of ∆µ a and ∆D was performed, which resulted in three volumetric images of changes in oxyhemoglobin, deoxyhemoglobin, and total hemoglobin (Eqs.19 and 20) in the human forehead.Each reconstruction consisted of 25 iterations of the conjugate gradient scheme and took 4 hours each on a Pentium III 550 MHz processor.All volumetric images pertaining to oxyhemoglobin, deoxyhemoglobin, and total hemoglobin changes were combined into three movies displaying temporal changes in oxyhemoglobin, deoxyhemoglobin, and total hemoglobin, respectively.These movies are shown in Figs.4-6.To assist in the visualization, three different views (coronal, parasagittal, and horizontal) are displayed.The different colors within the color scale represent isosurfaces.The isosurfaces for oxyhemoglobin, deoxyhemoglobin, and total hemoglobin concentration represent surfaces where a change in mM, relative to baseline, is observed.In Fig. 4, a movie is shown depicting the changes in deoxyhemoglobin during the Valsalva maneuver.As one views the time series an initial punctate rise in deoxyhemoglobin can be seen.These localized areas coalesce into two main zones where the level of deoxyhemoglobin continues to increase until a change of 0.008 mM is observed.Continuing, the zones begin to merge and a large area shows a 0.001 mM concentration increase.However, notice two inner areas with a 0.004 mM decreased deoxyhemoglobin concentration.Finally, towards the end of the Valsalva maneuver two zones are prominent.At the center of these zones a maximum change of 0.047 mM is reached.Note that the color scales do not increase by fixed step values.The chosen nonlinear scale better reflects the nonlinear changes in deoxyhemoglobin during the Valsalva maneuver.Also, since one isosurface may be imbedded within another, some isosurfaces appear to have changed color.In Fig. 5, a movie is shown depicting the distribution of oxyhemoglobin during the Valsalva maneuver.As one views the time series a diffuse 0.01 mM increase and localized 0.01 mM decrease in concentration from baseline is observed.Also, the decrease in oxyhemoglobin occurs deeper from the surface.As one sustains the Valsalva the increase in (C) 2001 OSA 10 September 2001 / Vol. 9, No. 6 / OPTICS EXPRESS 283 oxyhemoglobin continues while the decrease in oxyhemoglobin from baseline appears to fluctuate.Finally, two zones display a marked 0.12 mM increase in oxyhemoglobin concentration.
In Fig. 6, a movie is shown depicting the total volume change, as reflected by total hemoglobin, during the Valsalva maneuver.As one views the time series an initial 0.01 mM decrease in total hemoglobin can be seen approximately 0.5 cm below the surface.As the maneuver is continued a 0.02 mM increase is observed closer to the surface.Toward the end of the maneuver a maximum increase of 0.17 mM is reached.This maximum is located in a relatively small zone and represents the sum of the previous two images, even though the zone may be obstructed in the deoxyhemoglobin image.The present findings are consistent with previously reported effects concerning optical responses to hemodynamic perturbations.Specifically, those groups measuring functional tasks record optical signal changes of approximately 1-2%.The resulting changes in the concentration of oxyhemoglobin, as for example shown in the topographic maps by Franceschini et al. [4], range from -0.002 mM to +0.003 mΜ.In this study we have shown that the Valsalva maneuver leads to changes in the measured optical signal of up to 40%, which is about 20 to 40 times stronger that the changes observed by Franceschini et al.We find in our study that these stronger (up to 40 times) changes in optical signal reflect similarly stronger changes (up to 40 times) in oxyhemoglobin (+0.120 mM) as well.
Furthermore, we find that changes in oxyhemoglobin are stronger than changes in deoxyhemoglobin.This is consistent with findings by Hoshi et al. [40] and Franceschini et al.
[4].Hoshi et al. noted that oxyhemoglobin is the more sensitive parameter of changes in CBF, and the direction of changes in deoxyhemoglobin is determined by the degree of changes in venous blood oxygenation and volume.From the data that Franceschini et al. provides in [4] one can derive ratios of changes in oxyhemoglobin to deoxyhemoglobin in the range of 2 to 6.In our case, we find that the ratio of the observed maximum concentration change of oxyhemoglobin (0.13) to deoxyhemoglobin (0.05) is 2.6.

DISCUSSION
In order to interpret the physiological change seen in Fig. 4-6, it is necessary to understand the physiology of a Valsalva maneuver.During a Valsalva maneuver the subject performs a forced expiration against a closed glottis.This leads to an increase in intrathoracic pressure and via transmission on the vascular tree, an increase in arterial blood pressure and venous pressure [41].
The physiological effect of the Valsalva on both the peripheral and deep cortical vascular bed is well documented.For example, in some cases the Valsalva maneuver has been known to cause temporary cerebral ischemia and fainting [42], and has also been associated with aneurysm rupture and subarachnoid hemorrhage [43].More specifically, venous hypertension occurs over the valveless cerebral veins and an increased blood volume in the overlying area is observed.However, for periods of a few seconds, the initial phase of the Valsalva, the rise in venous pressure is transmitted proportionally to the cerebral spinal fluid (CSF), which diminishes transmural pressure reducing cerebral vascular resistance and maintaining blood flow.If intrathoracic pressure is held at positive values for long periods (remainder of Valsalva) CSF pressure starts to fall back towards normal levels and both arterial and venous pial vessels dilate.With continued elevation a cessation of pial blood flow and cortical blanching and dilation of sulcal vessels has been observed [44].
To date only two groups have used near-infrared light to study the hemodynamics during a Valsalva maneuver.Employing the NIRO 500 spectrometer (Hamamatsu Photonics, Japan) with an optode separation of 5-6 cm on the forehead, Matta et al. [45] concluded that both venous outflow and cerebral blood flow are affected by the Valsalva maneuver.This is in agreement with several studies using Transcranial Doppler [46][47][48] that measured middle cerebral artery velocity during the maneuver.On the other hand, in a recent study by Steinbrink et al. [49], in which the authors used a time resolved system, it was suggested that the Valsalva maneuver equally affects early and late photons.The authors attributed this effect to an absorption change in the upper layers of the head (i.e.scalp).However, as the authors also state, their analysis relies on a linear model that may be incapable of detecting deeper located phenomenon.
By using a MOBIIR scheme, the volumetric reconstruction, in principle, is capable of separating blood oxygenation effects from deoxygenation effects.Hence, one can suggest that the changes in oxygenation are associated with factors affecting cerebral blood flow, while the changes in deoxygenation may reflect venous pooling.Our findings seem to indicate that flow changes are the more dominant effect over the forehead during the Valsalva maneuver.Also, one might expect that the resulting cerebral blood flow changes would produce a more spatially uniform hemodynamic change, rather than the localized points of activation as observed initially in Fig. 4-6.However, these initial punctate changes may be explained by the interaction of light with both the superficial and deep vascular beds.Specifically, a superficial scalp vessel sitting close to a detector on the surface may mask the deeper phenomenon and produce images showing focal points of changing hemodynamics.As one maintains a sustained expiration against a closed glottis the deeper vascular changes contribute a larger portion to the recorded signal and deeper, more spatially uniform, images become visible.
Finally, it should be mentioned that brain motion might be partially responsible for the measured signal changes.The pressure changes caused by the Valsalva maneuver lead to an initial caudal displacement and subsequent cranial displacement of the neutral position of the brain.However, recent studies [50] have shown that most of the movement occurs at the level of the brain stem, with only minor movement occurring over the forehead.Interestingly, the displacement of the head may call into question the role played by CSF.Since, the CSF also plays a crucial role in maintaining adequate blood flow to the brain when venous pressure is altered with straining, or changes in posture, some of the optical signal changes may be due to the scattering changes that may be concomitantly occurring with the Valsalva maneuver.In future works, we plan to analyze the experimental data set using a transport based reconstruction algorithm, which will approach answering some of these questions.

SUMMARY
In this paper we reported on the first volumetric reconstructions of the optical properties in the human head.Using a 3-dimensional model-based iterative image reconstruction (MOBIIR) scheme we generated images of oxyhemoglobin, deoxyhemoglobin, and total hemoglobin from measurements on the human forehead during a Valsalva maneuver.
The MOBIIR scheme comprised of a forward model that predicts the detector readings based on an assumed spatial distribution of optical properties, an objective function Φ that compares predicted with measured signals, and a scheme that updates the assumed optical parameters for subsequent forward calculations.As a forward model for light propagation in the head, we used a three-dimensional, continuous wave diffusion equation with the Robin boundary conditions, which was solved numerically using a finite-element method.The objective function was defined to accommodate difference-measurement data, which compared measurements during the Valsalva maneuver to a rest state.As our updating scheme we employed an iterative conjugate gradient method that minimized the objective function.The necessary gradient of the objective function with respect to the optical properties was determined via the method of adjoint differentiation.Furthermore, to use the MOBIIR scheme in volumetric brain imaging a description of the external contour of the head (air-tissue boundary) and accurate knowledge of optode location is required.To this end, we have implemented a pragmatic photogrammetric technique.
The three-dimensional optical tomographic studies of the Valsalva maneuvers showed that a redistribution of hemoglobin does occur and that the strongest changes appear to be localized approximately 1.0 cm from the surface of the forehead.We found that changes in oxyhemoglobin are several times stronger that changes in deoxyhemoglobin.While further studies are needed to fully understand the hemodynamic processes that occur during a Valsalva maneuver, the protocol employed allowed us to determine the plausibility of using DOT to determine the global changes in oxyhemoglobin and deoxyhemoglobin in the human head.
In the future, volumetric MOBIIR-type techniques may help assess the dynamic cerebral autoregulatory mechanisms that maintain constant cerebral blood flow within a wide range of cerebral perfusion pressures.The clinical importance of assessing the cerebral autoregulatory mechanism lies in the protection of the brain from the sequelae of arterial hypotension and hypertension and ultimately in the diagnosis of cerebrovascular disease [35,47].
only once, independent on the number of sources.As was mentioned earlier, to get the , derivative is calculated, where (C) 2001 OSA 10 September 2001 / Vol. 9, No. 6 / OPTICS EXPRESS 277

( 21 )
This filter sets the value at a given time point x i to the average value of its n neighbors to the left and to the right.In this work we chose n = 2, which produced a smooth function without much distortion in the overall intensity profile.From Fig.2it can be seen that the Valsalva maneuver protocol is very reproducible from epoch to epoch.One can observe a strong drop in all detectors, and at the peak of the Valsalva maneuver the measurement intensities changed with respect to the rest period by up to 40%.Other source-detector combinations (not shown) lead to similar traces.

Fig. 2 :
Fig. 2: Time series of the experimental protocol at 760nm: 3 epochs -consisting of a Valsalva maneuver followed by a rest period (between dashed black lines).

Figure 3 :
Figure 3: Time series for first Valsalva maneuver epoch for source 3 (flat red trace) at 760nm.

Fig 6 .
Fig 6. (1.6MB) Time-series reconstruction showing change in blood volume as reflected by total hemoglobin in units of mM during the Valsalva maneuver; Coronal view (upper left), Parasagittal view (upper right), Horizontal view (lower left), scale (lower right).