1. Introduction

Diffuse optical imaging (DOI) is a non-invasive technology that uses low-levels of non-ionizing light in the range of 650-900nm to record changes in the optical absorption and scattering of tissue (reviewed in [1]). Optical imaging can be used to record human brain activity using a spatially arranged grid light emitters and sensors arranged on the surface of the scalp as shown in Fig. 1 . Each of these measurement combinations records light transmitted from the source position, through the head, and exiting beneath a detector position. Due to the high scattering properties of biological tissue, which result in a diffuse nature of migrating photons of light, each of these measurement pairs is sensitive to optical absorption changes within a diffuse volume of tissue lying between the source and detector position. Based on previous simulation studies, at a typical source-detector distance of around 3cm, diffuse optics can be used to sample down to about 3-5mm of the outer surface of the brain [2,3]. During functional activation studies, an optical probe like the one shown in Fig. 1 is positioned over the suspected region of brain activity. Stimulus evoked changes in blood flow to the brain result in changes in optical absorption due to oxy- and deoxy-hemoglobin, hereafter termed HbO₂ and Hb respectively. This results in a change in the signal in the measurements over this area of the brain, which can be converted to hemoglobin changes by the modified Beer-Lambert relationship [4] for each channel as seen in Fig. 1. As compared with functional MRI (fMRI), DOI is less costly, portable, and allows for a wider range of experimental scenarios because it does not require a dedicated scanner nor require the subject to lie supine. Moreover, optical imaging has the ability to resolve changes in both oxy- and deoxy-hemoglobin within the brain using multiple wavelengths of light, which can potentially lead to the ability to discriminate blood flow and oxygen metabolism changes [5]. Over the last several years, DOI has been used to explore brain activity in a variety of experimental conditions including infant development [6], mobility and balance [7], cognitive [8,9], and sensor-motor [10] systems to name only a few.

 

Fig. 1 Example of optical head cap and evoked response signals. An example of the optical imaging head cap is shown on the left. This cap contains light sources and detectors connected via fiber optics. Each of the nearest-neighbor measurement pairs records changes in the optical absorption of the tissue in between the two positions. During functional studies, the arrangement of sensors can be used to estimate the spatial distribution of oxy- (HbO₂), deoxy- (Hb), and total-hemoglobin (HbT) at each position. The image on the right shows an example data set showing activation in the frontal portion of the probe.

Download Full Size | PDF

One of the current challenges to optical imaging studies is the reconstruction of accurate images of the changes in brain activity. As reviewed in [1], optical imaging represents a generally underdetermined inverse problem, in which the image of interest has many more unknown parameters then measurements. This is also an ill-posed inverse problem, which means that a single unique solution cannot be obtained. Many methods based on regularization [1,11] and Bayesian models [12–14] have been proposed to improve the constraints on the optical inverse problem. To date however, these proposals have focused on improvement of images from a single imaging session or subject. Group-level analysis— the estimation of average brain activity image from a group of more than one subject or session, poses additional difficulties. Namely, since the optical inverse model for a single subject’s data is ill-posed, a solution needs to be selected that is consistent with both that single data set and with the images from other members of the group. Since the image reconstruction of any one subject’s data can yield a nearly infinite variety of admissible solutions to the inverse problem by virtue of its ill-posed nature, the reconstruction of individual solutions followed by averaging across subjects will compound the errors introduced by reconstruction and will lead to loss in group-level statistical effects due to inconsistencies between reconstructions.

To date, most of the group-level analysis in optical imaging is done either based on selection and averaging of regions-of-interest or based on direct averaging in the native measurement space (e.g. optical source-detector pairings) and requires careful reproducible placement of the optical probe. For example, some groups record optical positions relative to the international 10-20 system. The limitation of region-of-interest or per channel approaches is that it should be restricted to analysis within a specific group of subjects. This approach may lead to false results, for example, in the case of comparisons between two age groups or genders where brain anatomy may systematically differ. Region-of-interest analysis will also not suffice to capture differences in the location or extent of the brain activity. For example, due to the nature of the optical measurement model (forward model), an observed smaller change in optical signal could mean that there was either a lesser hemoglobin change or that the activation was deeper from the surface. In the same manner, the same optical signal could result if the location of the activity was more diffuse. Thus, region-of-interest based group analysis can be used to indicate that two groups are different, but cannot be used to elicit the exact nature of this difference— magnitude of change, location of signal, or anatomical differences between subjects (e.g. systematic errors in the placement of the optical probe).

In this work, we propose a random-effects model for simultaneously reconstructing optical images from multiple subjects to estimate group-level statistics. The key feature of this model is not to necessarily reconstruct a more accurate group-level image (although we will show that this does occur), but rather to directly test the null hypothesis that two groups of subjects can be modeled by the same image of brain activity. In other words, we can reject the hypothesis that two groups of subjects are identical if and only if no solution to the inverse problem can be found that matches all the data from both subjects. If a solution does exist, this solution is not unique as per the ill-posed nature of the inverse problem and therefore we still can’t uniquely say precisely how two groups are different in terms of a completely accurate reconstruction model. However, because forward model used in the image reconstruction accounts for both differences in subject anatomy and differences in the positioning of optical sensors between subjects, this model does allow us to (at least to the accuracy of our forward model) rule out that such differences are due to systematic anatomical or probe differences.

In our random-effects model, each subject is approximated as a perturbation (random-effect) from the group image. We will then concatenate all the subjects’ data into a larger inverse model to be solved. The image reconstruction model that we will use is based on the cortical surface approximation method, which we recently described in [15]. In this method, brain activity is estimated directly on the surface of the cortex. This surface is generated using the FreeSurfer program [16,17], which creates an inflated and anatomically registered representation of the surface of brain. This surface registration is a necessary requirement of our current model since it allows registration of the parameter (image) space between subjects. In comparison, volumetric registration as done in similar fMRI models would require masking of the superficial skin and cerebral spinal fluid (CSF) layers. In the context of our group reconstruction model, this could mean that a brain voxel from one subject might align in volume space with a CSF voxel in a second subject. Thus, the edges of the brain in the reconstruction would be distorted, which is where optical measurements are most sensitive to brain activity. In fMRI, this is not as detrimental of an issue as it is for optical data. Using the proposed surface based method, optical image reconstruction in the space of this inter-subject registered surface provides a means to perform group-level averaging while accounting for both differences in the anatomy of the brain and placement of the optical sensors between subjects. Subsequent to our random-effects model, we will also describe an adaptation of our original surface model [18] to model systemic physiology. In order to solve the group-level inverse problem, we describe a Bayesian inversion of this random-effects model, which uses Restricted Maximum Likelihood (ReML) [19] to optimize the use of prior information in the reconstruction model. Recently, we described the application of ReML methods to optical imaging [14].

In this work, we will present three numerical examples to demonstrate the utility of this model. In the first example, we will present the case of estimating functional activities of a single subject with data simulating the imperfect placement of the optical probe over multiple sessions. This example will allow us to examine how the group-analysis model is affected by the sensitivity profile of a nearest-neighbor optical probe while simplifying the problem by keeping the subject anatomy the same (within-subject averaging). In the second example, we consider the case of optical imaging across multiple subjects, where now in addition to variations in the optical probe sensitivity, additional errors are introduced by differences in brain anatomy from one subject to another. The final example will compare two subject populations where in addition to the issues of varying probe locations and brain anatomies, each group has a different average functional activity. Here, we will assess the ability of our model to detect differences in brain activity between two groups of subjects.

Theory

The theory and implementation of diffuse optical imaging is the subject of several recent reviews [1,20–22]. In this section, we will first briefly describe the optical forward model and review our cortical-surface model that we have previously described in [18]. We will then describe the implementation of our random-effect group model.

Optical forward model

Diffuse optical imaging is a non-invasive technology that measures changes in the absorption and scattering of tissue between a pair of light emitters and detectors placed on the surface of the scalp. At each of these emitter positions, light is sent into the tissue at two or more wavelengths. Due to the highly scattering nature of biological tissue, this light spreads as it enters the tissue. The propagation of light through tissue is often approximated by a diffusion-based random walk of the photons of light and can be modeled through Monte Carlo or finite elements methods. The ensemble path of the light depends on the structure of the underlying layers of the tissue, which for the head is the scalp, skull, cerebral spinal fluid, and gray/white matter brain tissue, which defines the volume sampled by each optical emitter and detector combination. During brain activity, regional changes in blood flow alter the concentration of blood and in turn change the absorption of the tissue. For the typically small changes in absorption associated with brain activity (and neglecting scattering changes that are assumed to be even smaller), the change in optical absorption (optical density; OD) at a given wavelength is approximated by a first-order linear expression given by

where A _i,j is the optical measurement (forward) model obtained from estimation of the diffusion of photons through the tissue and describes the summation of absorption values along the diffuse path traveled by the light going from a particular light emitter to a detector pair (i,j), HbO₂ and Hb are respectively the levels of oxy- and deoxy-hemoglobin, ε_HbX is the molar extinction coefficient for HbO₂ and Hb at a particular wavelength. For reasons that should be clear shortly, we have defined two random noise terms in this expression; υ is additive measurement space noise (e.g. instrument noise on a detector) and ω_X is an additive physiological noise in HbO₂ and Hb. For our model, the measurement noise (υ) will be assumed to be a zero-mean normally distributed random variable with independence between optical channels and wavelengths. The physiological noise (ω) will also be assumed zero-mean, but is not necessarily independent between HbO₂ and Hb. For a given optical probe arrangement which consists of several source-detector pairs, the set of linear equations corresponding to all pairs at all wavelengths can be concatenated into a single linear expression of the form

where the measurement contains all wavelengths and measurement pairs

and the parameters of interest model oxy- and deoxy-hemoglobin changes

Both the oxy- and deoxy-hemoglobin terms are actually vectors representing the change at each point (voxel) in the image. The forward operator H is created from the concatenation of all the individual source-detector measurement expressions and encapsulates both the forward model projection from the volume underneath the optical probe to the set of measurement channels and the modified Beer-Lambert relationship (conversion of hemoglobin changes to absorption changes at multiple wavelengths). Note that Eq. (2) is the typical optical forward model with the exception of the explicit distinction of the two noise terms.

Surface-based image reconstruction

In a recent publication, we described an approximation to the optical inverse problem, which modeled brain activity on the surface of the cortex using available structural magnetic resonance imaging (MRI) information [18]. In this approach, rather than modeling the changes in oxy- and deoxy-hemoglobin over the three-dimensional volume of the brain, we instead modeled the values at the two-dimensional surface boundary of the brain cortex. In that work, we numerically showed that this approach produced fairly accurate recovery of brain signals (generated from a proper volumetric simulation of activation) and we experimentally demonstrated that this approach produced images consistent with functional MRI in two separate experiments (visual stimulation and a frontal cortex task).

A schematic of the surface-reconstruction model is presented in Fig. 2 . In brief, structural MRI for a given subject was first segmented using the Massachusetts General Hospital’s FreeSurfer [16,17] and MNE software (www.nmr.mgh.harvard.edu) into skin, skull, cerebral spinal fluid, and gray/white brain regions. The FreeSurfer program was then used to perform a tessellation of the cortical surface and an expansion (non-linear morphing) of the surface to a spherical representation [17]. Important to this current work, this spherical surface is registered across subjects such that the same location on the sphere (surface coordinate) corresponds to the same anatomical area in each subject [17]. We described how we then constructed a basis set for image reconstruction using a set of spherical wavelets. This could be thought of as an extension of a typical wavelet representation of a two-dimensional image, but in our case that two-dimensional image is wrapped around a closed folded surface, which has no edges. In brief, this wavelet basis was based on the work by Schröder and Sweldens [23] which describes the construction of the spherical wavelets from multiple resolution mesh obtained via recursive subdivision of an icosahedron approximating a sphere. The mesh is created by repeatedly subdividing an initial mesh such that each resolves into four “child” triangles at each new subdivision level. This is accomplished by defining a midpoint at each edge, followed by connecting the new midpoint together, leading to a higher resolution mesh. With each additional wavelet level, the model is able to depict higher resolution details. The point-spread function for levels 3-6 of this cortical model is shown in Fig. 3 . The surface wavelet basis is used to reparameterize the image of oxy- and deoxy-hemoglobin at the surface of the brain cortex [18]. Substituting this change of basis into Eq. (2) gives

 

Fig. 2 Schematic of analysis route. The figure above schematically illustrates the steps in the construction of the surface-wavelet model. This model was previously described in [18]. Using a structural MRI, the brain is segmented and then inflated into a spherical representation using the Massachusetts General Hospital’s FreeSurfer tools (www.nmr.mgh.harvard.edu). The cortical surface of the registered brain is used to generate a set of two-dimensional wavelets that depict textures (e.g. brain activity images) on the surface. The segmented brain is also used to construct the optical forward model.

Download Full Size | PDF

 

Fig. 3 Multiple resolution surface meshes for head and brain. In this work, a two-layered nested model was used to depict signals originating in the brain and the scalp. Both layers were modeled using a spherical wavelet basis, which is constructed from repeated subdivision of an icosohedral mesh as described in [18]. In panel A, we show the mesh generated for the surface of the brain. The first column shows the mesh at four levels (3-6) of the model and the second column shows the blurring of a 2cm object on the surface of the brain. Panel B shows the mesh (levels 1-4) and smoothing effect of a 2cm object on the surface of the skin.

Download Full Size | PDF

where the subscript w denotes the parameters and noise in the wavelet domain. H_brain is a subset of the original forward model corresponding to only the columns containing the nodes on the surface of the brain.

The matrix operators W_brain and W_brain ⁻¹ are the forward (analysis) and inverse (synthesis) wavelet operators for the brain respectively. In the original work, we used a truncated wavelet basis, which acted like a spatial smoothing operator. In this work, we will be using the full (square) wavelet matrix and use the ReML method to impose spatial smoothing priors.

Model of systemic physiology

Our original cortical surface model only described brain activation. However, optical data is often (always) contaminated by systemic physiology from the scalp. To extend our original model, we propose to add a second layer corresponding to the surface of the scalp. In the same way that the surface of the brain was represented by the spherical wavelet basis, we can define a second surface corresponding to the extracted surface of the scalp. The purpose of this layer is to model the superficial noise in the optical signal due to systemic physiology. To avoid boundary effects to the forward model, we used the outer surface of the skull rather then the skin itself and used the MNE software to extract this boundary from the same anatomical MRI that was used for the cortical surface extraction and inflation. A schematic of the nested two-layer model is shown in Fig. 3.

In comparison to the brain surface, which was inflated and registered across subjects, the skin layer is subject-specific. The complete two-layered forward model in the wavelet domain for a single subject is given by

Note, that our model contains four sources of physiological noise (ω) corresponding to oxy- and deoxy-hemoglobin in the skin and brain layers in addition to a measurement noise (υ) term per channel and per wavelength of light measured. Since we have transformed the parameter space, the wavelet domain physiological noise is given by

Group-level random-effects model

The conventional approach to the optical inverse model would involve the inversion of Eq. (2) to create an image of brain activity for each subject. An averaging of the resulting images across the subjects or groups would typically follow the inverse of each individual subject’s model. This approach is currently used in analysis models for MEG and EEG [24]. The difficulty with this approach for optical data is because each of these inverse models is ill-posed (which means that there is no unique solution), the solution to one particular subject might be different from a second subject. However, this does not imply that the brain activity of the two subjects were actually different from each other. This leads to a high false discovery rate of inter-subject differences. Since the models are ill-posed, there are infinite equally valid solutions that explain the data. For this reason, group-level statistics for optical imaging are difficult and lead to a loss of statistical power when preformed this way. The goal of our proposed model was to create a group-analysis inverse model that attempts to find a solution to the larger inverse model that is consistent with all subjects’ data. In other words, our null hypothesis might be that two subjects (subject A and subject B) are not significantly different from each other. Thus, we will setup a single inverse model to address whether subject A and B can be modeled by the same inverse solution. Note, this does not require us to get the inverse solution correct, since the larger inverse model is still under-determined and ill-posed. Rather this states that if any solution exists that can explain both sets of data, then we cannot reject the null hypothesis. To do this, we will now describe a random-effects group analysis model for optical imaging.

In the random-effects model, the image for each subject is treated as a perturbation (random effect) from a group image. That is

where Δβ^Subject is assumed to be a zero-mean random variable. Our reconstruction model allows us to define a common parameter space across all the subjects by virtue of the anatomical registration of the MRI brains. Our surface brain model uses the FreeSurfer program to register the cortical surface of all the subjects and to create a mesh of the surface such that any node on the surface (entry of β) corresponds to the same anatomy across all the subjects. For example, the precentral gyrus is represented by the same entry in the parameter space for all subjects despite the fact that the physical dimensions, location, and folding of the brain may differ between subjects. We assume that structural anatomy and functional anatomy correspond to each other, which is also the basis of group analysis of functional MRI data in the FreeSurfer model. This registration is true of the brain parameters but not the skin parameters since the skin surface is generated differently and in general, the systemic physiology will differ between subjects. This feature of the model allows us to write Eq. (9).

The random-effects model is created by concatenating the measurement models from each subject and is given by the expression

where we have used the variable HW to denote the wavelet domain forward operator (HW = H·W⁻¹). In this example, Eq. (10) specifically describes the condition in which we have two subjects, however, this model can be easily generalized. The reconstruction of group average is given directly by the first set of entries in the parameter vector (${\beta }_{W,brain}^{Group}$) whereas the individual subjects are recovered by the sum of the group image and the specific perturbation term (e.g. ${\beta }_{W,brain}^{SubjectA}={\beta }_{W,brain}^{Group}+\Delta {\beta }_{W,brain}^{SubjectA}$). Differences in the two subjects can be tested via a t-test of the difference in the two perturbation terms (e.g. $p\left(\Delta {\beta }_{W,brain}^{SubjectA}-\Delta {\beta }_{W,brain}^{SubjectB}\ne 0\right)$.

In itself, the random-effects model does not help us to solve the inverse model. In fact, given the definition of the random-effect model in Eq. (10), we have succeeded in making the problem even more ill-posed since by adding the group term, we now have (N + 1) sets of parameters where N is the number of subjects. However, by placing the model in this framework we will soon demonstrate that we are able to make use of a Bayesian formulation and restricted maximum likelihood (ReML) methods as means to solve this model. Using the ReML model [14], we can estimate the covariance of the various noise terms. Specifically, the covariance of the random-effect is assumed to be zero-mean random variable. Through the ReML optimized weighted-least squares expression, non-zero values for the group parameters are more favored then non-zero random effects; meaning that when possible the model will load weight into the joint-group image over placing an equal weigh onto each of the individual subject perturbation images. For example in the case that two subjects are identical (${\beta }_{W,brain}^{SubjectA}={\beta }_{W,brain}^{SubjectB}$), loading of the group variable (${\beta }_{W,brain}^{Group}\ne 0$and $\Delta {\beta }_{W,brain}^{SubjectA}=\Delta {\beta }_{W,brain}^{SubjectB}=0$) is more favorable then placing loads on the random perturbations (e.g. ${\beta }_{W,brain}^{Group}=0$and $\Delta {\beta }_{W,brain}^{SubjectA}\ne 0$and $\Delta {\beta }_{W,brain}^{SubjectB}\ne 0$).

As written in Eq. (10), this model would describe a one-level statistical test in which the subject label (in this case “A” and “B”) is controlled for. By extension of this model, additional levels can be introduced, for example in the two level model consisting of four subjects (A, B, C, and D) with two subjects in each of two groups (I,II) a two-level model is constructed such that:

which would have a forward model of the form (ignoring skin models for simplicity):

Similarly, group-level fixed or mixed-effects models can be created by dropping off the columns corresponding to specific random-effects (effectively setting a specific$\Delta \beta =0$ at the subject or group level)

Restricted Maximum Likelihood (ReML) estimate of the inverse solution

In order to solve the group-level model, we will use a method called Restricted Maximum Likelihood (ReML) which is an approach used to optimize Bayesian inverse models. This approach has been described in numerous texts in the context of fMRI analysis (reviewed in [24]) and the inverse problem for magneto- and electro-encephalography (MEG/EEG) [25]. Recently, we described the application of this model to diffuse optical imaging [14]. In this section, we will briefly describe the concept of ReML and refer to previous work for more detailed mathematical treatment.

In a Bayesian framework, the weighted least-squares cost function is given by the expression

where $Y-\widehat{Y}$is the residual error of the model and $\beta -{\beta }_0$is the difference between the estimate of the model parameters and some a priori expectation of the parameters. We have used the notation ${\Vert X\Vert }^2{}_{C^{-1}}\equiv X^T{C}^{-1}X$, which is the weighted L2 norm. In lay terms, the Bayesian interpretation of Eq. (13) is that C_N should be the covariance of the measurement error (e.g. υ from Eq. (2) and C_P is the covariance of the intrinsic uncertainty in the parameters (e.g. ω from Eq. (2). Thus, Eq. (13) describes a weighted distance measure, which combines both actual measurements and prior expectations of the model. For readers familiar with a regularization interpretation, we note that the Bayesian model is a subset of the general class of linear quadratic regularization models (for example, Tikhonov regularization) in which the regularization penalty (e.g. $\lambda \cdot C$ in the term${\Vert X\Vert }^2{}_{\lambda \cdot C}$) is specifically interpreted to be the inverse of a covariance matrix under the Bayesian model

The solution to the linear model $Y=H\cdot \beta$for both the Bayesian model or the linear quadratic regularization model is given by the Gauss-Markov expression.

with

In the Bayesian interpretation, $C_{\beta |Y}$is the covariance of the parameters given the measurements and encapsulates both the intrinsic uncertainty in the parameters and the noise of the measurement process projected onto the parameter space.

In both the Bayesian and regularization interpretations, the Gauss-Markov expression requires an estimate of the terms C_N and C_P. Since these terms are specifically covariances in the Bayesian framework, approaches such as maximum likelihood can be used to optimize the values of these terms. It can be shown that the solution that provides the maximum likelihood of the data conditional on some set of hyperparameters (p(Y|λ)) can be found by maximizing the free energy expression(see [24] for derivation)

Note, that the free energy expression is similar to original weighted-least squares cost function expression (Eq. (13) with the addition of terms to penalize for large covariance terms.

The covariance of the model are parameterized such that

where the set of Λ’s are hyperparameters that will be estimated in the model and the Q terms are a set of covariance components that describe a priori expectations about the structure of the covariance model. We will later further detail these components for our specific model within the methods section of this paper. The method of ReML is to maximize Eq. (16) by solving for the parameters (β) and the hyperparameters (Λ). The approach used in the SPM software [24] and our current work is to use the expectation-maximization algorithm [26] in which the solution to Eq. (16) is obtained by iteratively solving for the parameters under the assumption of a fixed set of hyperparameters (via the Gauss-Markov expression) and then solving for the hyperparameters that maximize Eq. (16) under the assumption of a fixed set of parameters. This second step can be done by taking the derivative of free energy expression with respect to the vector of hyperparameters and setting this derivative equal to zero. Since this approach has been published several times, our purpose here is not to exhaustively detail this model. Refer appendix III of [24] for details of the derivation and implementation of this approach in the context of the fMRI and MEG models.

ReML solution to the random-effects DOI model

For our random-effects DOI model, we will use a hierarchical Bayesian formulation of the forward model. To summarize our previous sections in light of this, our optical model represents a three level model and is given by the set of equations

In this model, ${\beta }_{0|W,brain}^{Group}$is a prior on the expected value of the brain image for the group. For the purpose of this work, we will assume this prior is zero, which enforces a minimum norm prior on the group parameter space. However, throughout this text we will leave this term in the description since an obvious extension of this model would be to use a non-zero prior from previous NIRS or fMRI data. In the hierarchical model, there are four separate noise terms, which are all zero mean random variables; υ^Subject is the measurement noise, $\Delta {\beta }_{W,brain}^{Subject}$is the random-effect for each subject, ${\omega }_{W,skin}$is the physiological noise in the skin for each subject, and finally ${\omega }_{W,brain}$ is the uncertainty in the group image of brain activity. In a similar manner to the per subject random-effect model, a multiple groups comparison model can be created by adding a fourth level to the model and so on. We will define the following noise models

The free-energy expression for our specific model is given by

In order to solve this model, we will parameterize the covariance models such that

In the following methods section, we will describe the details of these covariance components.

2. Methods

In this section, we will describe the implementation of the group-analysis model, having presented the theory behind this model in the last section. We will also describe the methods used to generate the optical simulations used in this paper.

Implementation of ReML and Covariance models

ReML was used to solve the linear inverse model under multiple prior constraints. This model was implemented in Matlab (version 2010a) through a modification of the program spm_reml.m, which is part of the SPM-8 software program (http://www.fil.ion.ucl.ac.uk/spm/ [24]; ). In the situation of our cortical surface model, each subject’s image contains several thousand unknowns, but only several dozen measurements (our specific surface model contained ~7,500 parameters for the brain and ~1,000 parameters for the skin model for each of HbO₂ and Hb). In the case of the group-analysis model, there are a total of [(N + 1)*<number of brain parameters> + N*<number of skin parameters>] parameters where N is the number of subjects. For five subjects, this is greater the 200,000 parameters. Thus, this is a massively under-determined model whose dimensionality quickly exceeds the computational capacity of most computers. In order to solve this model numerically and computationally efficiently, we first reparameterize the model using a singular value decomposition (SVD). For the model $Y=H\cdot \beta$, a singular-value decomposition ($H=U\cdot S\cdot V^T$) can be used to define a new set of parameters such that

and the covariance of this new set of parameters is now given by

Importantly, in this new space, the size of the forward model (U) is square and of dimensions matching the number of measurements per subject times the number of subjects. Moreover, the hyperparameters of this space are identical to those of the original model. Thus, in practice, the SVD reparameterized model can be solved to yield the estimate of the hyperparameters in a very efficient manner using the much smaller matrix form. Once the hyperparameters have been determined, the original larger model can be solved directly in a single-step from the Gauss-Markov expression (Eq. (14). Since we had found that the ReML computation ran close to O(N^2) with the size of the model, by performing the SVD reparameterization we were able to cut the computation of the model from over several years of theoretical computational time to only minutes. The accuracy of this reparameterization was verified from both a theoretical and empirical prospective using smaller dimensional ill-posed DOI inverse models. This approach is similar to the source-space ReML models used for MEG reconstruction [25,27].

In order to constrain the covariance terms to have only positive variance (real-valued data), we reparametered the hyperparameters with an exponential function such that

The ReML code was iterated on the SVD-reparameterized model until the value of the (exponential) hyperparameter changed by less then 1E-5.

Covariance components

In our model, there are four types of noise/uncertainty terms (Eq. (21). Each of these noise terms was modeled as a linear combination of covariance components weighted by the hyperparameters. E.g.

Measurement noise (C_N). For each subject, two diagonal (identity) covariance components were used to model the measurement noise with one component corresponding to each of the two wavelengths of light used in the model (specifically 830nm and 690nm were simulated). Thus, the number of terms in the model of measurement noise was 2 x N where N was the number of subjects.

Subject-Parameter Uncertainty in the brain (C_B). For each subject, two diagonal (identity) covariance component terms were used to model oxy- and deoxy-hemoglobin changes respectively. The number of terms in the model of random effect was 2 x N where N was the number of subjects.

Subject-Parameter Uncertainty in the skin (C_S). An additional two diagonal covariance component terms were used to model oxy- and deoxy-hemoglobin changes in the superficial (skin) layer. These covariance components were used to bias the superficial image to lower spatial resolutions (spatial smoothing) by imposing that the variance of each level of the wavelet model (from lowest to highest spatial frequencies) was attenuated as $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\sigma }^{2\cdot (n-1)}$}\right.$ where n is the corresponding level of the wavelet model for those parameters in the model and σ is a smoothing factor. We used a value of σ = 15mm for this work. The number of terms in the model of skin physiology was 2 x N where N was the number of subjects.

Group-Parameter Uncertainty (C_G). Finally, the covariance of the group-level model was modeled by three covariance component terms. Two unity diaginal terms were used to model the oxy- and deoxy-hemoglobin images respectively. A third term with negative unity along the off-diagonal elements corresponding to covariance between oxy- and deoxy-hemoglobin was used to statistically model the expected anti-correlation between oxy- and deoxy-hemoglobin. Note that in the ReML model, all the covariance components are actually soft statistical (covariance) priors. In total there were 6xN + 3 covariance components for the model where N was the number of subjects.

Calculation of the optical forward model

The procedure for the calculation of the optical forward model is described in [18]. In brief, a structural MRI image (MPRAGE image collected on a Siemens 3T TRIO scanner or the Colin27 atlas [28] as noted) was segmented, cortically inflated and registered using the software tools in FreeSurfer [16,17]. The segmented head (skin, skull, cerebral spinal fluid, and gray/white brain matter) was used to construct a finite element model and the program NIRFAST [29] was used to compute the optical forward model at simulated wavelengths of 690nm and 830nm using the optical properties described in [18,30]. The two wavelength forward models were combined with the extinction coefficients of hemoglobin (Beer-Lambert) to generate the spectral forward model. The registered FreeSurfer surface (e.g. “sphere.reg”) was used to construct the spherical wavelet model as described in [18]. An additional wavelet model of the skin was generated from the outer surface of the skull, which was obtained using the MNE software extension to the Massachusetts General Hospital’s FreeSurfer program. We used the outer surface of the skull rather then the actual surface of the head because we found that the later suffered from artifacts due to the optical forward model right at the surface.

In order to test the group-level model, we preformed two types of simulations (later described as example I and II/III). In the first type of simulations, we used a single MRI brain (for this we used the Colin27 atlas [28]; ) and simulated the repositioning of an optical probe on the head with varied degrees of error (mean displacements of 0mm, +/−0.5mm, +/−1.5mm, +/−2.5mm, +/−5.0mm, +/−10.0mm, and +/−15.0mm [half the source-detector distance]). This set of simulations was intended to look at the errors introduced by the positioning of optical sensors. In the second set of simulations, we used MRI images from five separate subjects and positioned the optical probe on the head of each image according to the distance from the nasion, inion, and ears using an initial affine registration followed by relaxation onto the surface (see [31]). In this second set of simulations, we are able to look at both the effects of differing subject anatomy and positioning differences in the optical probe. In all simulations, we assumed that we knew the position of the probe on the head precisely, but that the probe positioning varied across subjects with random error.

In all simulations, we have used a nearest-neighbor probe as shown in Fig. 3. Although representing the majority of published optical studies, the nearest-neighbor probe is known to have areas of low sensitivity, particularly the volume directly beneath a source or detector position. With this type of probe, partial volume errors associated with the placement of the probe relative to the area of activation can be severe. Previous work by Joseph et al [32] provides evidence that nearest-neighbor geometries are inferior to high-density probe arrangements with overlapping (tomographic) measurement geometries. In our opinion, however, these high-density probes are not always practical depending on the area of interest, the hair color and thickness of the subject, and the capabilities of the instrumentation. Thus, in this work, we examined the performance of our group-level model in comparison to individual reconstructions for the case of the nearest-neighbor probe only. It is expected that the reconstruction of tomographic (e.g [32].) probe in the case of a single subject/session would be more robust to partial volume and probe positioning errors. Since the vast majority of DOI studies currently do not use tomographic probes, we feel that this study is currently justified.

Simulations

In order to demonstrate this model, we will go through three examples of increasing complexity. All simulations were preformed using a nearest-neighbor optical probe with a source-detector spacing of 32mm which we currently also use for experimental studies. The layout of the probe was shown in Fig. 1. Figure 4 shows the sensitivity of this probe to the underlying gyri and sulci in the brain. This probe has 4 source and 8 detector positions providing 15 measurement combinations at each of the two wavelengths (690nm and 830nm).

 

Fig. 4 Sensitivity of the optical probe. The sensitivity model for the optical probe used in this study (see Fig. 1) is shown panel A. This probe consisted of 4 source and 8 detector positions in a nearest-neighbor arrangement. The relative sensitivity map is shown on the inflated surface of the brain for the Colin27 atlas. The anatomical gyri and sulci locations are indicated in green and blue respectively in panel B. The sensitivity of the optical probe is near zero within the sulci of the brain.

Download Full Size | PDF

To simulate brain activity, a 7cm^2 homogenous target was generated in the center of the probe. The location of this activation was chosen for one of the subject’s brains as a reference and then placed in the registered inflated surface coordinate system. The location of the simulated activation was then transferred onto the individual brains of each simulated subject. A superficial model of systemic noise was also simulated. For each individual simulation, a noise image was created on the surface of the skin by randomly placing a 5cm radius inclusion on the skin layer followed by spatial low-pass filtering and renormalization to an amplitude of between 1 and 5μm for oxy-hemoglobin. Due to the low-pass filtering operation, the additive skin noise contained both positive and negative signed values. Deoxy-hemoglobin noise was added as −0.25 times the oxy-hemoglobin image. The noise in the skin layer was randomly positioned for each simulation. An example of a noise model used to represent the skin is shown in Fig. 5 .

 

Fig. 5 Skin noise model. In addition to brain activation, a low-spatial frequency noise image was added to the skin layer of the model. This was generated by placing an inclusion with a 5cm radius on the surface of the skin layer followed by spatial smoothing. The skin noise image was randomly placed for each simulation. An example noise image for oxy-hemoglobin is shown above.

Download Full Size | PDF

In the first example, we will show results from the simulation of a within-subject study in which we examine the effects of the repositioning of the DOI probe with random displacements (see Fig. 6 ). In this example, the anatomy and simulated brain activity were the same for each simulation but the location of the probe on the head, the additive systemic noise and the additive random measurement noise was different. In the second example, we simulated a multiple-subject study in which different anatomical volumes and repositioned DOI probe placements were used. For each subject, the same brain activation pattern was generated in the registered spherical space and then warped onto that subject’s individual cortical map (Fig. 7 ). Thus, each subject had slightly different location and extent of brain activity simulated depending on their actual anatomy. As before, additive systemic noise and the additive random measurement noise was different between subjects. In this example, we are examining the effects of both probe sensitivity and anatomical differences on the optical reconstruction. In the third and final example, we simulated a multiple-subject, multiple-session, and multiple-condition study. Here, five separate subjects were simulated in two sessions each. The simulated brain activity of the first and second sessions differed. In this example, we will look at the performance of the models to detect differences in brain activity subject to errors introduced by probe placement, subject anatomy, and systemic noise.

 

Fig. 6 Positions of the optical probe for within subject simulations. Optical data was simulated for displacements of the optical probe of 0.5-15mm. The probe was repositioned five times for each level of displacement error. In panel A, we show the procedure to register the optical probe on the head of a phantom using a Polhemus three-dimensional digitizer. Random perturbations were generated numerically from a registered template. In panel B, the placement of the five probes for the +/−15mm displacement simulations is shown. A different color marker represents each position. For simulations, the Colin27 [28] atlas was used.

Download Full Size | PDF

 

Fig. 7 Generation of inter-subject simulations. Optical data from multiple subjects was simulated by first generating an activation region on the registered atlas (spherical) map. This activation was then morphed back onto the individual cortical structures for five subjects via FreeSurfer. According to the registration of the cortical surface, the simulated activation in each of the five subjects represents the same anatomical area in the inflated space, however the extent and location in the three-dimensional (folded) brains varies between subjects.

Download Full Size | PDF

Comparison of models

We will compare three versions of the group-level model. First, images were individually reconstructed using a ReML model and then were averaged together in the image (parameter) domain. Overall statistical maps were constructed assuming independence between the individual subject images. In the second method, we used a mixed-effects model in which the random-effect term was left off allowing for direct reconstruction of the group average via joint-reconstruction. In the third method, we used the proposed random-effects model. For all models with the exception of those simulated with zero noise, a total of ten simulations were preformed with varied levels of random additive noise added to the measurement channels.

In all of our examples, we are performing a sort of inverse crime in that we are assumed to know the DOI probe position and its sensitivity model exactly. We have also assumed that the registration offered by FreeSurfer is exact in order to generate individual activation maps from the registered spherical atlas. In reality, we would expect errors in the determination of the position of the probe, errors in the registration of FreeSurfer, and inaccuracies of the optical forward model. The purpose of this current work is to explore whether a random-effects model can improve group-level analysis compared to individual subject level image reconstructions followed by averaging. Thus, the use of this simplified inverse model is justified to address this specific question. Future work to improve the optical forward model and optical registration will be needed to fully realize our model.

3. Results

Our proposed random-effects model out performed the approach of performing separate reconstructions for each subject followed by averaging in terms of both accuracy and statistical effect size. In this section, we will highlight three examples of this model going from the case of a within subject design to a between groups design. The later is more difficult since brain activation pattern, optical probe placement, and subject anatomy must be accounted for.

Example I. Within-subjects design

A known limitation of the DOI nearest-neighbor probe design is the location of regions of low sensitivity, in particular, directly beneath source or detector positions on the probe. Although tomographic probes using overlapping optical measurements provide more uniform sensitivity [32], the complicated setup and instrumentation of these probes is a limiting factor in their success. In general, most optical studies still use nearest-neighbor probe designs. In this first example, we simulated the repeated placement of an optical probe on the Colin27 atlas brain in order to examine the reconstruction errors introduced by a nearest-neighbor probe. We compared the individual reconstruction of each data set to joint reconstructions using a fixed-effect and random-effect model. A total of ten iterations of simulations were preformed for each of five levels of signal-to-noise ratios (infinite, 50:1, 20:1, 10:1, and 1:1) and each of seven mean displacement errors (0-15mm) of the probe. For each simulation, the probe was repositioned five times. A total of 1750 simulations were preformed.

In Fig. 8 , we show a single set of simulations for a probe displacement of +/15mm at a signal-to-noise level of 10:1. In Fig. 9 , the activation T-statistics for the same data is shown. As can be seen in both figures, compared to the individual reconstruction/average both the mixed and random-effects models produced much more localized activation estimates of oxy-hemoglobin that closely agreed with the simulated truth image. The effects size of the estimated activation was four-fold higher for both the fixed and random effects models. In the case of deoxy-hemoglobin reconstructions, all three methods produced comparable activation amplitudes, but the random and fixed effects models showed vast improvements to the statistical effect size. For this example, since the same signal-to-noise level using random measurement noise was simulated for all models, we don’t expect a large difference between the fixed- and random-effects analysis in this case. This does demonstrate that, as expected, the random-effects terms in the model are approaching zero since it is possible to model all five simulated data sets using a single group-level image.

 

Fig. 8 Intra-subject optical reconstructions from repeated measurements. The images above show the results of a simulation of intra-subject data with a random optical probe displacement of +/−15mm (also see Fig. 6). Data was simulated with a signal-to-noise ratio of 10:1. In column A, the simulated truth images are shown. The top/bottom rows show oxy- and deoxy-hemoglobin respectively. Column B shows the average of five individual reconstructions using a ReML model. Column C and D show the same data reconstructed using a fixed and random-effects model respectively.

Download Full Size | PDF

 

Fig. 9 Statistical effects sized from intra-subject optical reconstructions. The images above show the t-statistic for the data presented in Fig. 8. Based on Satterthwaite estimate [24] of the effective degrees-of-freedom for the model (1.8), the cutoff for significance of 0.05 is 3.1 standard deviations. Both the fixed and random effects images (columns C and D respectively) exceed this threshold for the area over the simulated region, but individual reconstruction and averaging results (column B) shows no significant activation.

Download Full Size | PDF

In Fig. 10 , we show bar graphs of the reconstruction accuracy and model fit to the data for the various reconstructions at each signal to noise level. At all five signal-to-noise levels the fixed and random-effects models were nearly identical. Both preformed better then the individual reconstruction followed by averaging method. The ReML method, which results in an empirically determined weighted least-squares model, adjusted the measurement noise term across the five signal-to-noise levels. As a result, the root mean squared measurement error ($\sqrt{{\Vert Y-\widehat{Y}\Vert }^2}$) approximately followed the additive level of measurement noise. In contrast, the model error ($\sqrt{{\Vert {\beta }_{True}-\beta \Vert }^2}$) was roughly the same across all signal-to-noise levels. This is because the parameter noise term (C_P) in the ReML model was approximately the same order of magnitude across all reconstructions. This is expected of the ReML model since differing levels of random measurement noise but not parameter noise was added in the simulations.

 

Fig. 10 Reconstruction errors in within subject simulation. In panel A, the residual error of the model is shown for each of the five signal-to-noise levels and for the individual, fixed, and random-effects model. In panel B, the model error (reconstructed group image minus the simulated image) for the brain is shown.

Download Full Size | PDF

Example II. Across-subjects design

In this second example, we simulated a between subjects study consisting of DOI signals simulated from five different subjects. Both the anatomy and location of the optical probe varied between subjects. Data was simulated based on the registered (inflated) cortical space. In Fig. 11 , we show the reconstructed maps of oxy-hemoglobin for the average of the five subjects and the group images for the fixed and random-effects models. In the multiple-subject model, the estimate of brain activity is directly obtained in the registered spherical space (inflated cortical surface). Thus, we have shown the group images on the surface of the registered surface (top row). We have also shown the image of this group estimate on the surface of one of the five subject’s brain. Similar to the intra-subject reconstruction case, the fixed and random-effects model were better localized compared to the independent reconstruction model. The deoxy-hemoglobin results (not shown) were consistent with those previously shown for the intra-subject example. In Fig. 12 , we show the T-statistics for the same set of data. As before, the fixed and random-effects models out preformed the individual reconstructions and the effects size of the random and fixed effects reconstructs was three to four times that of the individual reconstructions.

 

Fig. 11 Reconstructions of simulated data in inter-subjects study. The images above show the results of a simulation of inter-subject data (N = 5). Data was simulated with a signal-to-noise ratio of 10:1. In column A, the simulated truth images are shown. The top/bottom rows show the reconstructions of oxy-hemoglobin on the registered spherical surface and the group image projected onto one particular cortical surface. Column B shows the average of five individual reconstructions. Column C and D show the same data reconstructed using a fixed and random-effects model respectively.

Download Full Size | PDF

 

Fig. 12 Statistical effects of the inter-subjects study. The images above show the t-statistic for the data presented in Fig. 11. As before, the cutoff for significance of p<0.05 is 3.1 standard deviations. Again, both the fixed and random effects images (columns C and D respectively) exceed this threshold for the area over the simulated region, but individual reconstruction and averaging results (column B) shows no significant activation. Column A shows the simulated true image.

Download Full Size | PDF

In both the individual reconstructions and the random-effects model, an estimate of brain activity is obtained for both the group average and the individual subjects. In the case of the random-effects model, the individual’s image is obtained from Eq. (10). In Fig. 13 , we show the simulated images and reconstructions for the five individual subjects. The reconstruction errors for various signal-to-noise ratios are shown in Fig. 14 . In the case of the individual reconstructions, we can see that two of the subjects (A and E) had fairly good reconstructions while two of the other subjects (B and D) were underestimated. This reflects the uneven sensitivity of the nearest-neighbor probe and arose since the activation for these two subjects happened to be in an area of low sensitivity of the optical probe (refer to [32] for a more detailed treatment of this issue). In the case of the random-effect model, the data from all other subjects is being used to improve any one reconstruction. Thus, the reconstructions for subjects B and D are improved because additional information about the location of the activation was conveyed from the other three subjects in, which the probe was better positioned. As noted before, the issue of non-uniform sensitivity is particular to the nearest-neighbor probe.

 

Fig. 13 Reconstructed images for the five subjects. Both the individual reconstructions and the random effects model allow recovery of both the group averages (previously shown in Figs. 11 and 12) and the estimates for each individual subject. The individual subjects’ images are recovered from the random-effects model through Eq. (9) and are shown above for the same data presented in Fig. 10. The top row shows the location of the simulated activity on each subject’s brain. Note that the location and extent varies based on the individual folds of the brain. All five activation patterns were generated from the same registered cortically registered brain space (see Fig. 7). The middle row shows the results of the individual reconstructions and the bottom row shows the random-effects model.

Download Full Size | PDF

 

Fig. 14 Reconstruction errors in inter-subject simulation. In panel A, the residual error of the model is shown for each of the five signal-to-noise levels and for the individual, fixed, and random-effects model for the data presented in Figs. 11-13. In panel B, the model error (reconstructed group image minus the simulated image) for the brain is shown.

Download Full Size | PDF

Example III. Multiple-session across subjects design

In this final example, we examined the case of an across subjects design consisting of two simulated sessions which each had different brain activation. In this set of simulations, we examined the performance of the models to detect differences in activation from one session to the other. We compared reconstructions using a three-level and four-level random-effects model. The three-level model was identical to the random-effects model described in the previous two sections (refer to Eqs. (18-20) for the definitions of the levels). This model consisted of a group term and five subject perturbation terms. Using this model, the two groups of subjects were examined separately and thus two independent random-effects models were solved. In the four-level model (see Eq. (12), there is an overall term, two group-level perturbation terms, and five subject-level perturbation terms. The two groups were simulated with amplitudes of 8 and 4μM with a small displacement in the location of the activation areas. The results of the image reconstructions are shown in Fig. 15 . The t-statistics for these reconstructions are shown in Fig. 16 . Both versions of the model preformed nearly equally for the estimation of the amplitude of the images. Both amplitudes for the two groups were accurately estimated. Although the amplitudes for these two images were correct and the difference between the images was around 4μM, we were unable to detect the difference in spatial location of the two groups. This was not too surprising since this remains an ill-posed inverse model. The t-statistics for the individual group images are comparable for both reconstructions as well. We were able to detect significant amplitude differences between the two groups, but we unable to detect the small difference in the position of the activity.

 

Fig. 15 Image reconstruction of between groups study. The images above show the results of a simulation of two groups of subjects (N = 5/5). The simulated activation in group I and II was 8μm and 4μm respectively for oxy-hemoglobin (column A) and −2μm and −1μm for deoxy-hemoglobin (not shown). Data was simulated with a signal-to-noise ratio of 10:1. Columns B and C show the reconstructions of the two groups using a three and four-level random-effects model respectively. In the three-level model (B), the two groups were reconstructed independently. In the four-level model (C), both groups (10 subjects total) were reconstructed simultaneously.

Download Full Size | PDF

 

Fig. 16 Statistical effects of the inter-groups study. The images above show the t-statistic for the data presented in Fig. 15. The cutoff for significance of p<0.05 is 2.2 standard deviations (Satterthwaite correction). For both group's the estimated activation was significant for both the three- and four-level random effects models (columns B and C respectively). The difference image was only significant in the case of the four-level random-effects model (Column C). Column A shows the simulated true image.

Download Full Size | PDF

4. Discussion

In this work, we have described a random-effects model that uses joint reconstruction of multiple sessions of optical data to directly test the null hypothesis that two groups of subjects are identical. We are able to reject this hypothesis if at least one solution to the joint reconstruction model can be found that models both sets of data. Because our model also uses data from all of the subjects to estimate group and individual (rand-effects) images, we were able to demonstrate that the joint group-level model offers improvement of the reconstruction accuracy of both group and individual level estimates. In particular, we presented a series of simulated examples to show that this model improves the ability to localize brain activity compared to the approach of reconstructing each subject individually. Our model also improved the recovered effects size for the reconstruction by three to four-fold over the averaging of individual reconstructions. In many cases, this improvement marked the difference between being able to statistically reject the null hypothesis (e.g. no significant activation detected). In the random-effects model, both the group image and the individual results are obtained whereas in the fixed-effect model only the group average is retained. Thus, the random-effects model still allows single-subject information, but in effect uses all the other subjects to improve the individual estimates.

Comparison of reconstruction models

Our data supports an advantage of a multi-subject joint image reconstruction compared to individual reconstructions followed by averaging. However, in most cases we found that the fixed and random-effects models were very similar for many of the simulations. We failed to find any situations in which the random-effects model was significantly better then the fixed effects one. We believe that the differences between the two models may become more apparent under more situations of the model, which should be examined in future work. In particular, if the physiological and measurement noise across subjects or groups is random and of comparable magnitudes, we would expect both of these models to perform equally. If instead, one or few of the subjects had been outliers (either in terms of brain activation or physiological noise levels), the random-effects model would be expected to perform better. The advantage of the random-effects model over the fixed effect model is that the random-effects model allows recovery of both the group and individual level data whereas only the fixed effects model only returns the group average result. Whether or not the individual images are of interest depends on the nature of the study and we not that the random-effects model has substantially more unknown parameters.

Model assumptions and limitations

We have made three major assumptions in this model. First, we have assumed brain activation to be constrained to the surface of the cortex. Our previous work [18] examined the justification and limits of this assumption. In short, we had found that cortical reconstructions closely matched both simulated volumetric data and experimental data collected with concurrent fMRI-DOI. In this work, we simulated activation directly on the surface of the brain and thus committed an inverse crime of sorts. In the original work, we had examined the performance of cortical reconstructions from volumetric simulations and had found that the cortical surface constraint resulted in slight under-estimation of the magnitude of the activation. The basis of the FreeSurfer multiple-subject registration is the alignment of the cortical surfaces rather then the brain volumes. Because of this, in the context of this work, it is impossible to generate registered volumetric simulations. In other words, we cannot simulate a single volumetric image of brain activation for the group that can be morphed into the individual brain volumes. Thus, the limitation of our cortical surface assumption depends on the accuracy of the cortical approximation.

A second assumption of our model was that the position of the DOI probe could be determined and registered to the MRI structural image with certainty. In our experimental models, we use a Polhemus three-dimensional scanner (as shown in Fig. 2) and an initial affine registration based on landmarks on the head followed by an automated relaxation algorithm to move the probe positions to the surface of the head whilst keeping the distances between the sensors along the surface fixed. Based on initial work, we believe that this registration is accurate to around 2-3mm, but we do not know to what extent an error in the probe positioning of this size will have on the reconstruction accuracy of the DOI data.

Finally, in this work we have modeled superficial physiology using a low-pass filter bank of wavelets, which existed on a surface around the skull. This approach is conceptually similar to work done in Huppert et al [13] where DOI and fMRI brain data was modeled by a low-resolution voxel basis set in the superficial layers. Here, we have modeled brain and skin changes using two concentric surfaces derived from the anatomy.

5. Conclusions

We have demonstrated how the registration and simultaneous reconstruction of optical data from multiple subjects can be used to improve group-level comparisons of functional data. In particular, our approach allows us to account for variations in the optical probe placement and underlying structural anatomy, which could confound analysis across subjects or experimental sessions. Because our method is based on the reconstruction of optical data within a common registered space (cortical surface space), the group-level reconstructions allow for the estimate of images that are simultaneously consistent with all of the subjects' data. We showed that this approach was more robust to partial volume and probe positioning errors and offered improvements to the effects size of the estimated group-level images.