*i*) the spatial resolution is <200 μm for NA ≤0.2 or focal plane depth ≤300 μm. (

*ii*) More than 97% of the signal comes from the top 500 μm of the tissue. (

*iii*) For activated columns with lateral size larger than spatial resolution, changing numerical aperature (NA) and focal plane depth does not affect depth sensitivity. (

*iv*) For either smaller columns or large columns covered by surface vessels, increasing NA and/or focal plane depth may improve depth sensitivity at deeper layers. Our results provide valuable guidance for the optimization of optical imaging systems and data interpretation.

## 1.

## Introduction

In two-dimensional (2-D) camera-based optical brain imaging either on an exposed neo-cortex or through a thinned skull, a collimated beam (often in the visible spectrum) incident onto the brain surface propagates inside the tissue. Depending on the optical properties being imaged (absorption or fluorescence), the optical imaging system collects either diffusely reflected or fluorescent light from the brain, forming 2-D images. Here, we disregard the effects of scattering changes because they are negligible in the visible spectrum. Both absorption [such as with intrinsic optical imaging (IOI), which images the absorption of the hemoglobin] and fluorescence imaging [such as with voltage-sensitive dye imaging (VSDI)] have been used widely to map neo-cortical functions *in vivo*.^{1, 2, 3, 4, 5, 6} For example, local blood oxygenation and total hemoglobin concentration often change with neuronal activity in the cortical tissue, leading to a change in the absorption of incident light, and therefore, a difference in diffusely reflected light.^{3, 5, 6} Hence, images formed by reflected light can reveal cortical function. On the other hand, fluorescence imaging may reveal different aspects of neuronal activity depending on the dyes employed. For example, voltage-sensitive dyes reflect the electric activities of neurons. Both approaches have been used extensively in studying the functional maps and spatiotemporal dynamics of the neuronal activity of the somatosensory, visual, and auditory cortexes, as well as the olfactory bulbs in animal models such as mice, rats, cats, and monkeys.^{3, 4, 5, 6, 7, 8}

Tissues at different depths of the cortex contribute to the image on the camera [i.e., the measured signal (the intensity of each pixel of the image) is a weighted sum of the real signal from these tissues]. We use the term “depth sensitivity” to describe this weighting function. Because of light scattering in the cortex, even a single point in the brain tissue will not result in an image consisting of a diffraction-limited spot but a much larger region. The size of this region is determined not only by the configuration of the optical imaging system [numerical aperture (NA) and focal plane depth] but also by the optical properties of the tissue (such as the scattering and absorption coefficients). We use the term “spatial resolution” to characterize the size of this region due to an array of points throughout the cortical tissue. In general, both the spatial resolution and depth sensitivity may depend on the optical configuration of the imaging system. Hence, in order to optimize the imaging system and better interpret the experimental data and compare the results among different functional imaging modalities (including IOI, fluorescence imaging, and nonoptical-based approaches, such as functional magnetic resonance imaging), we need to characterize the spatial resolution and depth sensitivity as a function of NA and focal plane depth. The light scattering and absorption in biological tissues make it difficult to do so; however, Polimeni. have estimated the spatial resolution of IOI for a specific optical configuration using Monte Carlo simulation with tissue parameters describing the mammalian cortex.^{9} Their result of 240 μm is consistent with the experimental result of Orbach and Cohen.^{10} Other authors have mentioned a much smaller spatial resolution, though without offering evidence to support it.^{3} In summary, spatial resolution has been reported for only a few specific optical configurations and there is a lack of information for depth sensitivity.

Propagation of light in a highly scattering medium such as the cortex is complicated and therefore, difficult to model by directly solving Maxwell equations. On the other hand, both absorption and fluorescence imaging, where the polarization and interference effects can be ignored, can be rigorously modeled with the radiative transfer equation (RTE). The analytical solution of RTE is known for only a limited number of simple experimental configurations, however, and it fails when we need to study light propagation near the boundaries or the light sources, as is the case for many 2-D optical brain imaging experiments. Monte Carlo simulation traces the trajectories of individual photons and offers a flexible yet rigorous approach to photon transport in highly scattering tissues. Since Wilson and Adam first introduced Monte Carlo simulation into the field of laser tissue interactions, it has been widely used to simulate light transport in tissues.^{11, 12, 13, 14}

We have employed Monte Carlo simulation to calculate the spatial resolution and depth sensitivity of both the absorption- and fluorescence-based 2-D optical imaging under optical configurations typically encountered in camera-based optical brain imaging on the neo-cortex. In Sec.2, we present physical models of the problem and describe our method to reduce the number of Monte Carlo simulations to improve the efficiency of our calculation. Here, we have employed the RTE under the first Born approximation. This enables us to perform the Monte Carlo simulations using homogeneous tissue configuration and thus greatly reduce computational time. In Sec. 3, we show numerical results for spatial resolution and depth sensitivity under various combinations of NA and focal plane depth. In Sec. 4, we discuss the application of our results in optical brain-imaging experiments *in vivo*, especially, how to optimize NA and focal plane depth differently to satisfy different requirements of spatial resolution, depth sensitivity, and signal level.

## 2.

## Methods

Figure 1 illustrates the 2-D camera-based optical imaging system. A collimated light beam perpendicularly illuminates the brain tissue. For absorption imaging, diffusedly reflected light (marked by arrows) is collected by lens L_{1} and imaged onto a camera by lens L_{2}.^{15} When there is no functional activation, the brain tissue scatters and absorbs light with homogeneous refractive index *n*
_{0}, the anisotropic factor *g*
_{0}, scattering, and absorption coefficients, μ_{s0} and μ_{a0}, respectively; thus, the image on the detector plane is uniform. When a certain cortical region is functionally activated, its local optical properties, such as the absorption coefficient, will often be altered. Such spatial disturbance of the otherwise homogeneous optical property will lead to a small perturbation in the image on the detector, resulting in a nonuniform image. Therefore, in a typical functional imaging experiment, images are recorded both during and before the functional activation and the difference image is obtained as δ*I*. For fluorescence imaging, fluorescence is imaged, instead of reflected light, by inserting a bandpass filter (BPF) in front of the detector. Thus, the difference image would reflect changes in local optical properties, such as the absorption coefficient and/or the fluorescent quantum efficiency η_{0}. Ideally, this difference image would reflect the local variation of the optical properties at a particular depth inside the tissue, which, in turn, would represent the functional activation pattern of the tissue. As stated in the Introduction, due to both the scattering nature of cortical tissue and the limitations of the imaging system, this difference image is a “blurred” representation of the spatial disturbance of the local optical properties throughout the cortical depth. In fact, in 2-D optical imaging experiments, we do not have enough information to reconstruct the complete 3-D variation of the local optical properties from the difference images. Nevertheless, these difference images offer valuable information as to the spatial extent and magnitude of the functional activation provided that we understand the spatial resolution and depth sensitivity of the 2-D optical imaging method, which we will define and calculate for both absorption and fluorescence imaging.

Figure 1 illustrates a “4*f* configuration” where the two lenses have the same focal length *f* and the same distances between the consecutive planes. To be precise the distances between the front focal plane of L_{1}, L_{1}, the back focal plane of L_{1} (the same as the front focal plane of L_{2}), L_{2}, the back focal plane of L_{2} (the same as the detector plane) are all *f*. Hence, the front focal plane of L_{1} in an optically clear medium is imaged directly onto the detector with a magnification *M* = 1. Because of its simplicity and wide application, we will mainly consider the 4*f* configuration in this paper.^{15} The potential deviation of our results from those of different configurations will be addressed briefly in Sec. 4. In many functional imaging experiments, the optical system can be varied by either changing the NA of L_{1} or placing different depths of the cortical tissue at the front focal plane of L_{1}, as shown in Fig. 1. We will discuss the influence of these changes on spatial resolution and depth sensitivity.

In order to define the spatial resolution for absorption imaging in the cortex, which often has columnar organization, we consider a single array of point absorbers along the optical axis and evenly distributed throughout the cortical depth [Fig. 2a]. Each point absorber acts as a local perturbation and results in a difference image with circular symmetry [Fig. 3a]. The summation of these images is then detected by the camera [Fig. 3b]. The spatial resolution can be defined as the full width at half maximum (FWHM) of that profile [Fig. 3c]. As displayed here, the spatial resolution indicates the minimal separation between two neighboring functionally activated columns that can be spatially resolved by the imaging system. To define the depth sensitivity, we consider multiple arrays of point absorbers evenly distributed throughout the cortical depth and with a cross-sectional area of *L*
^{2}, where *L* is the lateral size of the functionally activated region inside the tissue [Fig. 2b]. These arrays of point absorbers act as the local perturbation. Because the point absorbers at all depths contribute to the difference image, the total intensity variation at one location (or pixel) on the detector δ*I*(*x*
_{d}, *y*
_{d}) is the summation of all the individual intensity variation, δ*I*
_{z}(*x*
_{d}, *y*
_{d}), where *z* indicates the depth [Figs. 4a and 4b]. Therefore, we define depth sensitivity as the ratio of the contribution from an individual depth *z* to the total contribution from all depths, i.e., δ*I*
_{z}/δ*I* [Fig. 4c]. Similarly, we can define spatial resolution and depth sensitivity for fluorescence imaging by replacing the absorbers with fluorescent targets inside the tissue. Both spatial resolution and depth sensitivity play important roles in accurately interpreting the data and in comparing data obtained using different imaging modalities and across species.

As can be seen, we need to know the corresponding difference images for absorption and fluorescence imaging, respectively, in order to calculate the spatial resolution and depth sensitivity. This requires an understanding of light propagation in a highly scattering medium, such as the cortex. As discussed earlier, the RTE provides an accurate description of this problem. However, there is no analytical solution to the transport equation in the situation of 2-D optical imaging, where we are interested in light propagating from near to within a few scattering lengths from the surface. Fortunately, in a typical functional imaging experiment, the local perturbation to the optical properties is relatively weak: only a few percent of that of the background. Hence, we can simplify the problem by using the first Born approximation to the RTE followed by Monte Carlo simulations. Below we will discuss the cases for absorption imaging and for fluorescence imaging.

## 2.1.

### Absorption imaging

The total absorption coefficient can be written as a sum of a homogeneous background component μ_{a0} and a spatially varying perturbation δμ_{a0}(*r*), thus, μ_{a}(*r*) = μ_{a0} + δμ_{a}(*r*). As discussed above, we use the first Born approximation to the RTE and express the perturbed radiance
[TeX:]
$\delta I(r_{\rm d}, \hat \Omega _{\rm d})$
$\delta I({r}_{\mathrm{d}},{\widehat{\Omega}}_{\mathrm{d}})$
, the intensity flowing in direction
[TeX:]
$\hat \Omega _{\rm d}$
${\widehat{\Omega}}_{\mathrm{d}}$
and arriving at position *r*
_{d} on the detector plane, as follows:

## Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta I(r_{\rm d}, \hat \Omega _{\rm d}) = - \int\!\!\int\delta\mu_{\rm a}(r)I_0(r,\hat \Omega)G(r_{\rm d}, \hat \Omega _{\rm d} |r,\hat \Omega)d\hat \Omega d^3 r, \end{equation}\end{document} $$\delta I({r}_{\mathrm{d}},{\widehat{\Omega}}_{\mathrm{d}})=-\int \int \delta {\mu}_{\mathrm{a}}\left(r\right){I}_{0}(r,\widehat{\Omega})G({r}_{\mathrm{d}},{\widehat{\Omega}}_{\mathrm{d}}|r,\widehat{\Omega})d\widehat{\Omega}{d}^{3}r,$$*r*and

*r*

_{d}are the positions inside the tissue and on the detector plane, respectively, and [TeX:] $\hat \Omega$ $\widehat{\Omega}$ and [TeX:] $\hat \Omega _{\rm d}$ ${\widehat{\Omega}}_{\mathrm{d}}$ are the unit vectors in the directions that the light travels inside the tissue and toward the detector plane, respectively.

^{16, 17, 18, 19}[TeX:] $I_0 (r,\hat \Omega)$ ${I}_{0}(r,\widehat{\Omega})$ is the radiance inside the homogenous tissue described by optical properties of μ

_{s0}, μ

_{a0}, n

_{0}, and

*g*

_{0}due to collimated light illumination. [TeX:] $G(r_{\rm d}, \hat \Omega _{\rm d} |r,\hat \Omega)$ $G({r}_{\mathrm{d}},{\widehat{\Omega}}_{\mathrm{d}}|r,\widehat{\Omega})$ is the Green function solution for the transfer equation for light emitting from a point

*r*inside the same homogeneous tissue in the direction of [TeX:] $\hat \Omega$ $\widehat{\Omega}$ to a point

*r*

_{d}on the detector in the direction of [TeX:] $\hat \Omega _{\rm d}$ ${\widehat{\Omega}}_{\mathrm{d}}$ . Next, we apply Eq. 1 to the corresponding δμ

_{a}(

*r*) and obtain the spatial resolution and depth sensitivity, respectively. Note that the intensity at a particular location

*I*(

*r*) is the radiance arriving at that location integrated over all the solid angles, i.e.,

## Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} I(r) = \int {I(} r,\hat \Omega)d\hat \Omega . \end{equation}\end{document} $$I\left(r\right)=\int I(r,\widehat{\Omega})d\widehat{\Omega}.$$*r*

_{d}, δ

*I*(

*r*

_{d}) can be written as

## Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta I(r_{\rm d}) = - \int\!\!\int\delta\mu_{\rm a}(r)I_0(r,\hat \Omega)G(r_{\rm d} |r,\hat \Omega)d\hat \Omega d^3 r. \end{equation}\end{document} $$\delta I\left({r}_{\mathrm{d}}\right)=-\int \int \delta {\mu}_{\mathrm{a}}\left(r\right){I}_{0}(r,\widehat{\Omega})G\left({r}_{\mathrm{d}}\right|r,\widehat{\Omega})d\widehat{\Omega}{d}^{3}r.$$*r*inside the homogeneous tissue in the direction of [TeX:] $\hat \Omega$ $\widehat{\Omega}$ to a point

*r*

_{d}on the detector and is related to [TeX:] $G(r_{\rm d}, \hat \Omega _{\rm d} |r,\hat \Omega)$ $G({r}_{\mathrm{d}},{\widehat{\Omega}}_{\mathrm{d}}|r,\widehat{\Omega})$ by [TeX:] $G(r_{\rm d} |r,\hat \Omega)\break = \int_0^{2\pi } G (r_{\rm d}, \hat \Omega _{\rm d} |r,\hat \Omega)d\hat \Omega _{\rm d}$ $G\left({r}_{\mathrm{d}}\right|r,\widehat{\Omega})={\int}_{0}^{2\pi}G({r}_{\mathrm{d}},{\widehat{\Omega}}_{\mathrm{d}}|r,\widehat{\Omega})d{\widehat{\Omega}}_{\mathrm{d}}$ .

## 2.1.1.

#### Spatial resolution

To calculate the spatial resolution, we consider a single array of point absorbers evenly distributed throughout the cortical depth and with an absorption coefficient of μ_{a0} + Δμ_{a0}, where Δμ_{a0} is the extra absorption due to the functional activation; thus,

## Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta \mu _{\rm a} (r) = \sum\limits_{0 \le z{}_i \le 1\,{\rm mm}} {\Delta \mu _{{\rm a}0} \delta (x,y,z - z_i)}, \end{equation}\end{document} $$\delta {\mu}_{\mathrm{a}}\left(r\right)=\sum _{0\le z{}_{i}\le 1\phantom{\rule{0.16em}{0ex}}\mathrm{mm}}\Delta {\mu}_{\mathrm{a}0}\delta (x,y,z-{z}_{i}),$$*x*,

*y*,

*z*–

*z*

_{i}) is the Dirac δ function representing a point absorber at (0, 0,

*z*

_{i}) inside the tissue. In our calculation, each voxel is 10 μm in size, and we consider depths from 0 to 1 mm below the surface because the light intensity reaching the detector from tissues deeper than this is negligible for the wavelength range (in the visible up to ~670 nm) typically used in both absorption and fluorescence imaging. Therefore, there are 100 identical point absorbers [Fig. 2a]. Note that [TeX:] $I_0 (r,\hat \Omega)$ ${I}_{0}(r,\widehat{\Omega})$ in Eq. 1 only depends on depth

*z*due to the collimated illumination and the homogeneous background optical properties; thus, [TeX:] $I_0 (r,\hat \Omega) = I_0 (z,\hat \Omega)$ ${I}_{0}(r,\widehat{\Omega})={I}_{0}(z,\widehat{\Omega})$ . Combining Eqs. 4, 3, we obtain

## Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta I(r_{\rm d}) = - \Delta \mu _{{\rm a0}} \sum\limits_{0 \le z \le 1\,\,{\mathop{\rm mm}\nolimits} } {\int {I_0 (z,\hat \Omega)G(r_{\rm d} |0,0,z,\hat \Omega)d\hat \Omega } } . \end{equation}\end{document} $$\delta I\left({r}_{\mathrm{d}}\right)=-\Delta {\mu}_{\mathrm{a}0}\sum _{0\le z\le 1\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\mathrm{mm}}\int {I}_{0}(z,\widehat{\Omega})G\left({r}_{\mathrm{d}}\right|0,0,z,\widehat{\Omega})d\widehat{\Omega}.$$*I*(

*r*

_{d}) is the difference image caused by the single array of point absorbers. [TeX:] $G(r_{\rm d} |0,0,z,\hat \Omega)$ $G\left({r}_{\mathrm{d}}\right|0,0,z,\widehat{\Omega})$ is the Green function solution for the transport equation for light emitting from a point (0, 0,

*z*) inside the homogeneous tissue (where a point absorber is located) in the direction of [TeX:] $\hat \Omega$ $\widehat{\Omega}$ to a point

*r*

_{d}on the detector. Therefore, we can interpret δ

*I*(

*r*

_{d}) as follows [Fig. 5a]: collimated light illuminates the homogeneous tissue and leads to a radiance of [TeX:] $I_0 (z,\hat \Omega)$ ${I}_{0}(z,\widehat{\Omega})$ arriving in the direction of [TeX:] $\hat \Omega$ $\widehat{\Omega}$ at the point absorber located at (0, 0,

*z*). The point absorber then acts as a “negative light source” with strength proportional to [TeX:] $ - \Delta \mu _{{\rm a}0} I_0 (z,\hat \Omega)$ $-\Delta {\mu}_{\mathrm{a}0}{I}_{0}(z,\widehat{\Omega})$ that propagates back toward the detector. We then obtain the intensity variation at the detector due to this point absorber alone by integrating the radiance over 4π [i.e., [TeX:] $ - \Delta \mu _{{\rm a}0} \int {I_0 (z,\hat \Omega)G(r_{\rm d} |0,0,z,\hat \Omega)d\hat \Omega }$ $-\Delta {\mu}_{\mathrm{a}0}\int {I}_{0}(z,\widehat{\Omega})G\left({r}_{\mathrm{d}}\right|0,0,z,\widehat{\Omega})d\widehat{\Omega}$ ]. Hence, the total intensity variation at the detector caused by all the point absorbers is simply the summation of each individual contribution, as described by Eq. 5.

We employ Monte Carlo simulation to calculate both
[TeX:]
$I_0 (z,\hat \Omega)$
${I}_{0}(z,\widehat{\Omega})$
and
[TeX:]
$G(r_{\rm d} |0,0,z,\hat \Omega)$
$G\left({r}_{\mathrm{d}}\right|0,0,z,\widehat{\Omega})$
using the homogenous optical properties μ_{s0}, μ_{a0}, n_{0}, and *g*
_{0}. We need one Monte Carlo simulation for
[TeX:]
$I_0 (z,\hat \Omega)$
${I}_{0}(z,\widehat{\Omega})$
, however, 100 different simulations for
[TeX:]
$G(r_d |0,0,z,\hat \Omega)$
$G\left({r}_{d}\right|0,0,z,\widehat{\Omega})$
, one for each point absorber located at a different depth *z*. Because Monte Carlo simulation is time consuming, it would be advantageous to reduce the number of calculations. To this end, we take advantage of the translational invariance of the tissue, the principle of reciprocity, and the 4*f* configuration to further simplify Eq. 5 into

## Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \hspace*{-2pc}\delta I(x_{\rm d}, y_{\rm d}, z_{\rm d}) &=& - \Delta \mu _{a0} \sum\limits_{0 \le z \le 1\,{\rm mm}}\nonumber \\ && \hspace*{-.0pc}\times {\int\!\! {I_0 (z,\hat \Omega)G(x_{\rm d}, y_{\rm d}, z, - \hat \Omega |0,0,z_{\rm d})d\hat \Omega } } \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill \delta I({x}_{\mathrm{d}},{y}_{\mathrm{d}},{z}_{\mathrm{d}})& =& -\Delta {\mu}_{a0}\sum _{0\le z\le 1\phantom{\rule{0.16em}{0ex}}\mathrm{mm}}\hfill \\ & & \phantom{\rule{0.0pt}{0ex}}\times \int {I}_{0}(z,\widehat{\Omega})G({x}_{\mathrm{d}},{y}_{\mathrm{d}},z,-\widehat{\Omega}|0,0,{z}_{\mathrm{d}})d\widehat{\Omega}\hfill \end{array}$$*x*,

*y*,

*z*) and traveling in a direction [TeX:] $ - \hat \Omega$ $-\widehat{\Omega}$ . It results from a point source located at the center of the detector plane (0, 0,

*z*

_{d}) and thus can be obtained from only one Monte Carlo simulation. Equation 6 shows that the difference image can be obtained from only two Monte Carlo simulations: one for [TeX:] $I_0 (z,\hat \Omega)$ ${I}_{0}(z,\widehat{\Omega})$ and the other for [TeX:] $G(x,y,z, - \hat \Omega |0,0,z_{\rm d})$ $G(x,y,z,-\widehat{\Omega}|0,0,{z}_{\mathrm{d}})$ .

## 2.1.2.

#### Depth sensitivity

Here, we consider multiple arrays of point absorbers evenly distributed throughout the cortical depth (0–1 mm) and with a cross-sectional area of *L*
^{2} [Fig. 2b], thus,

## Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta \mu _{\rm a} (r) = \sum\limits_{\scriptstyle - L/2 \le x_i \le L/2 \hfill \atop {\scriptstyle - L/2 \le y_j \le L/2 \hfill \atop \scriptstyle 0 \le z_k \le 1\,\,{\rm mm} \hfill}} \Delta \mu _{{\rm a}0} \delta (x - x_i, y - y_j, z - z_k). \end{equation}\end{document} $$\delta {\mu}_{\mathrm{a}}\left(r\right)=\sum _{\genfrac{}{}{0pt}{}{{\scriptstyle -L/2\le {x}_{i}\le L/2}}{\genfrac{}{}{0pt}{}{{\scriptstyle -L/2\le {y}_{j}\le L/2}}{{\scriptstyle 0\le {z}_{k}\le 1\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\mathrm{mm}}}}}\Delta {\mu}_{\mathrm{a}0}\delta (x-{x}_{i},y-{y}_{j},z-{z}_{k}).$$*z*

_{d}). In addition, we take advantage of the principle of reciprocity and obtain the expression for intensity variation due to a single activated layer centered at depth

*z*

_{0}and with a thickness of Δ

*z*, δ

*I*

_{z}

_{0}(0, 0,

*z*

_{d}), and the whole functionally activated column δ

*I*(0, 0,

*z*

_{d}). Therefore,

## Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \delta I_{z0} (0,0,z_{\rm d}) &=& - \Delta \mu _{{\rm a}0} \sum\limits_{\scriptstyle - L/2 \le x \le L/2 \hfill \atop {\scriptstyle - L/2 \le y \le L/2 \hfill \atop \scriptstyle z_0 - \Delta z/2 \le z \le z_0 + \Delta z/2 \hfill}}\nonumber\\ && \times {\int {I_0 (z,\hat \Omega)G(x,y,z,\hat \Omega |0,0,z_{\rm d})d\hat \Omega } } \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill \delta {I}_{z0}(0,0,{z}_{\mathrm{d}})& =& -\Delta {\mu}_{\mathrm{a}0}\sum _{\genfrac{}{}{0pt}{}{{\scriptstyle -L/2\le x\le L/2}}{\genfrac{}{}{0pt}{}{{\scriptstyle -L/2\le y\le L/2}}{{\scriptstyle {z}_{0}-\Delta z/2\le z\le {z}_{0}+\Delta z/2}}}}\hfill \\ & & \times \int {I}_{0}(z,\widehat{\Omega})G(x,y,z,\widehat{\Omega}|0,0,{z}_{\mathrm{d}})d\widehat{\Omega}\hfill \end{array}$$## Eq. 9

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta I(0,0,z_{\rm d}) = \sum\limits_{0 \le z_0 \le 1\,\,{\rm mm}} {\delta I_{z0} (0,0,z_{\rm d})} . \end{equation}\end{document} $$\delta I(0,0,{z}_{\mathrm{d}})=\sum _{0\le {z}_{0}\le 1\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\mathrm{mm}}\delta {I}_{z0}(0,0,{z}_{\mathrm{d}}).$$*I*

_{z}

_{0}(0, 0,

*z*

_{d})/δ

*I*(0, 0,

*z*

_{d}).

## 2.2.

### Fluorescence imaging

In fluorescence imaging, a photon is first absorbed by a fluorescent target and subsequently emitted at a longer wavelength with certain fluorescent quantum efficiency. The functional activation can cause a local perturbation of the absorption coefficient, the fluorescence efficiency, or a combination of the two. The net result is a change of the local fluorescent intensity. For simplicity, we treat the tissue as a homogeneous scattering medium with optical properties of μ_{s0}, μ_{a0}, n_{0}, and *g*
_{0} with its fluorescent quantum efficiency varying between η_{0} and η_{0} + δη(*r*) corresponding to without and with the functional activation, respectively. This local variation of the fluorescent quantum efficiency will lead to a difference fluorescent image δ*I*. Similar to that described in the case of absorption imaging, we can find this difference image as follows [Fig. 5b]: collimated light propagates inside the tissue and arrives at position *r*, where a fluorescent target is located with an intensity of *I*
_{0}(*r*). The fluorescent target then emits incremental fluorescence due to the functional activation with an intensity δη(*r*)*I*
_{0}(*r*). This extra fluorescence propagates back toward the detector as described by the Green function *G*(*r*
_{d}|*r*), resulting in an incremental fluorescent intensity at position *r*
_{d} of the detector, δη(*r*)*I*
_{0}(*r*)*G*(*r*
_{d}|*r*). The total incremental fluorescent intensity at position *r*
_{d}, δ*I*(*r*
_{d}), is the summation of the individual contributions from all the fluorescent targets inside the tissue

## Eq. 10

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta I(r_{\rm d}) = \int {\delta \eta (r)I_0 (r)G(r_{\rm d} |r)d^3 r} . \end{equation}\end{document} $$\delta I\left({r}_{\mathrm{d}}\right)=\int \delta \eta \left(r\right){I}_{0}\left(r\right)G\left({r}_{\mathrm{d}}\right|r){d}^{3}r.$$*r*,

*I*

_{0}(

*r*), and the Green function

*G*(

*r*

_{d}|

*r*) are related to the radiance [TeX:] $I_0 (r,\hat \Omega)$ ${I}_{0}(r,\widehat{\Omega})$ and Green function [TeX:] $G(r_{\rm d}, \hat \Omega _{\rm d} |r,\hat \Omega)$ $G({r}_{\mathrm{d}},{\widehat{\Omega}}_{\mathrm{d}}|r,\widehat{\Omega})$ as follows:

## Eq. 11

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} I_0 (r) = \int_0^{4\pi } {I_0 (r,\hat \Omega)d} \hat \Omega \end{equation}\end{document} $${I}_{0}\left(r\right)={\int}_{0}^{4\pi}{I}_{0}(r,\widehat{\Omega})d\widehat{\Omega}$$## Eq. 12

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} G(r_{\rm d} |r) = \int_0^{4\pi } {\int_0^{2\pi } {G(r_{\rm d}, \hat \Omega _{\rm d} |r,\hat \Omega)d\hat \Omega _{\rm d} d\hat \Omega } }. \end{equation}\end{document} $$G\left({r}_{\mathrm{d}}\right|r)={\int}_{0}^{4\pi}{\int}_{0}^{2\pi}G({r}_{\mathrm{d}},{\widehat{\Omega}}_{\mathrm{d}}|r,\widehat{\Omega})d{\widehat{\Omega}}_{\mathrm{d}}d\widehat{\Omega}.$$Comparing Eqs. 3, 10, we can see that the difference image of the absorption imaging is sensitive to the angular distributions of both the radiance [TeX:] $I_0 (r,\hat \Omega)$ ${I}_{0}(r,\widehat{\Omega})$ and the Green function [TeX:] $G(r_{\rm d}, \hat \Omega _{\rm d} |r,\hat \Omega)$ $G({r}_{\mathrm{d}},{\widehat{\Omega}}_{\mathrm{d}}|r,\widehat{\Omega})$ , whereas that for fluorescence imaging only depends on the total intensity. The negative sign in Eq. 3 means that the increase in tissue absorption reduces the intensity of the image on the detector plane, leading to a negative difference image, while the positive sign in Eq. 10 means that an increase in the fluorescence quantum efficiency increases the intensity of the image on the detector plane, leading to a positive difference image.

## 2.2.1.

#### Spatial resolution

Here, we consider a single array of fluorescent targets distributed from 0 to 1 mm below the surface with a depth profile described by *A*(*z*); thus,

## Eq. 13

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta \eta (r) = \sum\limits_{0 \le z_i \le 1\,\,{\rm mm}} {\Delta \eta _0 A(z)\delta (x,y,z - z_i)} . \end{equation}\end{document} $$\delta \eta \left(r\right)=\sum _{0\le {z}_{i}\le 1\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\mathrm{mm}}\Delta {\eta}_{0}A\left(z\right)\delta (x,y,z-{z}_{i}).$$*z*direction, depending on the staining process.

^{20, 21, 22}Hence, we choose to perform our calculation for two types of situations, one with

*A*(

*z*) being uniform [Fig. 2c] and the other with

*A*(

*z*) being nonuniform with the maximal dye concentration appearing 200 μm below the surface [Fig. 2d], as reported by Ferezou. in their VSDI experiments.

^{20}

Similar to absorption imaging, *I*
_{0}(*r*) only depends on depth *z*, thus, *I*
_{0}(*r*) = *I*
_{0}(*z*), and we can simplify Eq. 10 into

## Eq. 14

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \delta I(x_{\rm d}, y{}_{\rm d},z_{\rm d}) = \sum\limits_{0 \le z \le 1\,\,{\rm mm}} {\Delta \eta _0 A(z)I_0 (z)G(x_{\rm d}, y_{\rm d}, z|0,0,z_{\rm d})} .\nonumber\\ \end{eqnarray}\end{document} $$\begin{array}{c}\hfill \delta I({x}_{\mathrm{d}},y{}_{\mathrm{d}},{z}_{\mathrm{d}})=\sum _{0\le z\le 1\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\mathrm{mm}}\Delta {\eta}_{0}A\left(z\right){I}_{0}\left(z\right)G({x}_{\mathrm{d}},{y}_{\mathrm{d}},z|0,0,{z}_{\mathrm{d}}).\end{array}$$*G*(

*x*

_{d},

*y*

_{d},

*z*|0, 0,

*z*

_{d}) is the light intensity at position (

*x*

_{d},

*y*

_{d},

*z*) inside the tissue originating from a point source at the center of the detector (0, 0,

*z*

_{d}). It can be obtained as the radiance [TeX:] $G(x_d, y_d, z, - \hat \Omega |0, 0, z_d)$ $G({x}_{d},{y}_{d},z,-\widehat{\Omega}|0,0,{z}_{d})$ integrated over a solid angle of 4π as follows:

## Eq. 15

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} G(x_{\rm d}, y_{\rm d}, z|0,0,z_{\rm d}) = \int_0^{4\pi } {G(x_{\rm d}, y_{\rm d}, z, - \hat \Omega |0,0,z_{\rm d})d} \hat \Omega . \end{equation}\end{document} $$G({x}_{\mathrm{d}},{y}_{\mathrm{d}},z|0,0,{z}_{\mathrm{d}})={\int}_{0}^{4\pi}G({x}_{\mathrm{d}},{y}_{\mathrm{d}},z,-\widehat{\Omega}|0,0,{z}_{\mathrm{d}})d\widehat{\Omega}.$$## 2.2.2.

#### Depth sensitivity

Here, we consider multiple arrays of fluorescent targets evenly distributed throughout the cortical depth (0–1 mm) and with a cross sectional area of *L*
^{2} [Fig. 2b]; thus,

## Eq. 16

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \delta \eta (r) = \sum\limits_{\scriptstyle - L/2 \le x_i \le L/2 \hfill \atop {\scriptstyle - L/2 \le y_j \le L/2 \hfill \atop \scriptstyle 0 \le z_k \le 1\,\,{\rm mm} \hfill}} \Delta \eta _0 A(z)\delta (x - x_i, y - y_j, z - z_k)\nonumber\\ \end{eqnarray}\end{document} $$\begin{array}{c}\hfill \delta \eta \left(r\right)=\sum _{\genfrac{}{}{0pt}{}{{\scriptstyle -L/2\le {x}_{i}\le L/2}}{\genfrac{}{}{0pt}{}{{\scriptstyle -L/2\le {y}_{j}\le L/2}}{{\scriptstyle 0\le {z}_{k}\le 1\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\mathrm{mm}}}}}\Delta {\eta}_{0}A\left(z\right)\delta (x-{x}_{i},y-{y}_{j},z-{z}_{k})\end{array}$$*z*

_{d}). We obtain the expression for the intensity variations due to a single activated layer centered at depth

*z*

_{0}and with a thickness of Δ

*z*and the whole functionally activated column as δ

*I*

_{z}

_{0}(0, 0,

*z*

_{d}) and δ

*I*(0, 0,

*z*

_{d}), respectively. Hence,

## Eq. 17

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta I_{z0} (0,0,z_{\rm d}) = \Delta \eta _0 \sum\limits_{\scriptstyle - L/2 \le x \le L/2 \hfill \atop {\scriptstyle - L/2 \le y \le L/2 \hfill \atop \scriptstyle z_0 - \Delta z/2 \le z \le z_0 + \Delta z/2 \hfill}}\!\!\!\!\!\!\!\!\! {A(z)I_0 (z)G(x,y,z|0,0,z_{\rm d})} \end{equation}\end{document} $$\delta {I}_{z0}(0,0,{z}_{\mathrm{d}})=\Delta {\eta}_{0}\sum _{\genfrac{}{}{0pt}{}{{\scriptstyle -L/2\le x\le L/2}}{\genfrac{}{}{0pt}{}{{\scriptstyle -L/2\le y\le L/2}}{{\scriptstyle {z}_{0}-\Delta z/2\le z\le {z}_{0}+\Delta z/2}}}}A\left(z\right){I}_{0}\left(z\right)G(x,y,z|0,0,{z}_{\mathrm{d}})$$## Eq. 18

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta I(0,0,z_{\rm d}) = \sum\limits_{0 \le z_0 \le 1\,\,{\rm mm}} {\delta I_{z0} (0,0,z_{\rm d})} . \end{equation}\end{document} $$\delta I(0,0,{z}_{\mathrm{d}})=\sum _{0\le {z}_{0}\le 1\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\mathrm{mm}}\delta {I}_{z0}(0,0,{z}_{\mathrm{d}}).$$*I*

_{z}

_{0}(0, 0,

*z*

_{d})/δ

*I*(0, 0,

*z*

_{d}).

We calculate
[TeX:]
$I_0 (z,\hat \Omega)$
${I}_{0}(z,\widehat{\Omega})$
and
[TeX:]
$G(x,y,z,\hat \Omega |0,0,z_{\rm d})$
$G(x,y,z,\widehat{\Omega}|0,0,{z}_{\mathrm{d}})$
using a publicly available Monte Carlo code, “tMCimg,” developed by Stott and Boas (see Refs. 19 and 23). It allows arbitrary initial position and direction of the photon and records the radiance within each voxel of the tissue, which are necessary for our calculations. The tissue cross-sectional area, depth, and voxel size are 2 × 2 mm^{2}, 5 mm, and 10 μm, respectively. To calculate
[TeX:]
$I_0 (z,\hat \Omega)$
${I}_{0}(z,\widehat{\Omega})$
, the initial direction of the photon is perpendicular to the tissue surface while the initial position is determined by two independent random numbers so that it is uniformly distributed within 2 × 2 mm^{2}. To calculate
[TeX:]
$G(x,y,z,\hat \Omega |0,0,z_{\rm d})$
$G(x,y,z,\widehat{\Omega}|0,0,{z}_{\mathrm{d}})$
, photons from a point source at the detector strike the tissue surface at a zenith angle θ [from 0 deg (perpendicular to the surface) to sin^{–1} (NA)] and an azimuth angle ψ (from 0 deg to 2π), both of which are in turn determined by two independent random numbers. The initial positions are determined by θ, ψ, and focal plane depth. The main advantage of tMCimg over the widely used code “MCML” developed by Wang and Jacques (see Ref. 14) lies in its ability to perform simulations in a medium with arbitrary boundaries and spatial variation in the optical properties. It takes a Xeon 2.26-GHz CPU 2 h to compute 10^{8} photons propagating in a volume of 2 × 2 × 5 mm^{3} with a voxel size of 10 × 10 × 10 μm^{3}.

In particular, We choose μ_{s0}= 35 mm^{–1}, μ_{a0} = 0.27 mm^{–1}, n_{0} = 1.4*,* and *g*
_{0} = 0.9, which are close to those that describe the gray matter of the mammalian cortex at 633 nm and can represent the tissue properties in the visible range up to 670 nm.^{9, 24} The radiance inside the tissue due to the collimated illumination
[TeX:]
$I_0 (z,\hat \Omega)$
${I}_{0}(z,\widehat{\Omega})$
does not depend on the NA and focal plane depth of the imaging system; hence, we need only one Monte Carlo simulation for it. The radiance inside the tissue due to the point source at the center of the detector
[TeX:]
$G(x,y,z,\hat \Omega |0, 0, z_{\rm d})$
$G(x,y,z,\widehat{\Omega}|0,0,{z}_{\mathrm{d}})$
depends on the NA and focal plane depth of the imaging system therefore, we must calculate
[TeX:]
$G(x,y,z,\hat \Omega |0,0,z_{\rm d})$
$G(x,y,z,\widehat{\Omega}|0,0,{z}_{\mathrm{d}})$
for each optical configuration. Here, we choose typical ranges of NA (0.1, 0.2, and 0.4) and focal plane depth (20, 100, 200, 300, 400, 500, and 600 μm) used in functional imaging experiments; therefore, 21 calculations will be needed for
[TeX:]
$G(x,y,z,\hat \Omega |0,0,z_{\rm d})$
$G(x,y,z,\widehat{\Omega}|0,0,{z}_{\mathrm{d}})$
. From the radiances
[TeX:]
$I_0 (z,\hat \Omega)$
${I}_{0}(z,\widehat{\Omega})$
and
[TeX:]
$G(x,\, y,\, z,\,\hat \Omega |0,\, 0,\, z_{\rm d})$
$G(x,\phantom{\rule{0.16em}{0ex}}y,\phantom{\rule{0.16em}{0ex}}z,\phantom{\rule{0.16em}{0ex}}\widehat{\Omega}|0,\phantom{\rule{0.16em}{0ex}}0,\phantom{\rule{0.16em}{0ex}}{z}_{\mathrm{d}})$
, we then obtain the intensities *I*
_{0}(*z*) and *G*(*x*
_{d}, *y*
_{d}, *z*|0, 0, *z*
_{d}) using Eqs. 11, 15. Subsequently, we get the corresponding difference images for both absorption and fluorescence imaging using Eqs. 6, 8, 9, 14, 17, 18. Finally, we calculate the spatial resolution and depth sensitivity using the corresponding difference images. The results are presented in Sec. 3.

## 3.

## Results

## 3.1.

### Spatial profile of the difference Image Widens and the Spatial Rresolution Worsens with the Increase of NA and Focal Plane Depth

Figure 6 depicts the spatial profiles of the difference images due to the presence of a single array of point absorbers (or fluorescent targets) for both absorption and fluorescence imaging and under different optical configurations. In particular, Fig. 6a shows the spatial profiles corresponding to NA = 0.2 and focal plane depths at 100 (blue color), 300 (red color), and 500 μm (black color), respectively, for absorption imaging. As the imaging system focuses deeper into the tissue, the spatial profile broadens. Figure 6b shows the spatial profiles corresponding to three different NAs [0.1 (blue), 0.2 (red), and 0.4 (black)] and focal plane depth = 300 μm for absorption imaging. As the imaging system collects more light that exits the tissue with large angles, the spatial profile broadens. Figures 6c and 6(d) depict the corresponding profiles for fluorescence imaging, which exhibit the same trend. We explain the results of Fig. 6 as follows. Assume we have calculated the difference image on the tissue surface generated by the array of point absorbers of Fig. 2a. We would expect it to have features similar to that of Fig. 3b, thus a similar profile along the radial direction depicted in Fig. 3c. Now, let us examine how this profile would be imaged onto the detector plane by the 4*f* system. If the focal plane depth is 0, then a point on the tissue surface will be imaged onto a point on the detector plane (assuming there is no aberration). Therefore, the detector plane would have a similar profile. However, if the focal plane depth is *d*, then an area with a radius ~*d* tan[sin^{–1}(NA/*n*
_{0})] will contribute to a point on the detector; thus, the profile on the detector would be “blurred.” The larger the focal plane depth *d* and NA are, the more blurring it will induce, thus leading to a wider profile and a lower spatial resolution (as can be seen in Fig. 7).

The uniform [straight line, Figs. 6c and 6d] and nonuniform [dashed line, Figs. 6c and 6d] distributions of the fluorescent targets have similar profiles, with the exception that the profile of the nonuniform distribution has slightly longer tails than that of the uniform distribution under the same optical configuration. As discussed below, only fluorescent targets within the top ~500 μm of the tissue have a major effect on the spatial profile. Figures 2c and 2d illustrate that, within the top 500 μm, more fluorescent targets lie in the deeper layers in the nonuniform distribution, resulting in more light scattering and thus longer tails in this case.

The spatial profiles of the absorption imaging have longer tails than their fluorescence imaging counterparts under the same optical configuration. Equation 5 indicates that radiance
[TeX:]
$I_0 (z,\hat \Omega)$
${I}_{0}(z,\widehat{\Omega})$
from the incidence light arriving at the point absorber at (0, 0, *z*) and flowing in the direction
[TeX:]
$\hat \Omega$
$\widehat{\Omega}$
. The point absorber absorbs it and subsequently acts as a negative light source with an emission that propagates toward the detector in the same direction
[TeX:]
$\hat \Omega$
$\widehat{\Omega}$
[Fig. 5a]. Equation 14 describes a similar process for fluorescence imaging except that after the radiance
[TeX:]
$I_0 (z,\hat \Omega)$
${I}_{0}(z,\widehat{\Omega})$
is absorbed, the fluorescent target emits isotropic fluorescence [Fig. 5b]. Because we are interested in light propagating from near to within a few scattering lengths from the surface and the scattering is highly forward (anisotropic factor *g* is 0.9), the point absorbers act as anisotropic sources that emit radiation preferentially toward deeper tissue layers. Consequently, this light would travel deeper into the tissue and thus experience more scattering as compared to the emission from fluorescent targets that act as isotropic light sources. Therefore, the spatial profile in the former case has longer tails.

Figure 7 shows the spatial resolution of absorption (a) and fluorescence imaging (b) versus focal plane depth under different NAs. In both cases, the spatial resolution worsens with the increase of the NA and focal plane depth. Note that the spatial resolutions of absorption and fluorescence imaging are similar, with the exception that those of absorption imaging are slightly worse under the same optical configurations, which is consistent with their spatial profiles. In general, as we vary the NA from 0.1 to 0.4 and the focal plane depth from 20 to 300 μm, we change the spatial resolution from <20 to 200 μm, which is still sufficient to resolve many individual functional columns. For example, the horizontal size of many vertebrate columns ranges from 50 to 800 μm.^{25} We can compare our result to the spatial resolution obtained by Polimeni using the following tissue parameters: μ_{s0}= 35.4 and 53.2 mm^{−1}, μ_{a0}= 0.27 and 0.22 mm^{–1}, n_{0} = 1.4*,* and *g*
_{0} = 0.94 and 0.82 for gray and white matter, respectively; and NA ~0.4 and focal plane depth = 300 μm.^{9} They used two lenses coupled directly instead of the 4*f* configuration. They obtained a spatial resolution of 240 μm for the summed contributions within a tissue depth of 200–500 μm.^{9} This is similar to the 200 μm, which we obtained for NA = 0.4, focal plane depth = 300 μm, and a tissue depth of 200–500 μm. The difference can be attributed to the slight variation in the tissue parameters and the difference between their lenses being directly coupled versus our 4*f* configuration.

The profiles shown in Fig. 6 have larger tails than that of a Gaussian with the same FWHM, which is consistent with that reported by Polimeni.^{9} Furthermore, Figs.
6c, 6d, 7b show that modifying the density distribution of the fluorescent dyes across the cortical depth from Figs. 2c to 2d can result in a spatial profile with longer tails and yet a slightly smaller FWHM width (i.e., the spatial resolution). This suggests that, in addition to spatial resolution, it is also important to consider a cross talk between signals generated in adjacent cortical columns, which is larger in a case of diffusely reflected light than in the case of Gaussian-like signal profiles.

## 3.2.

### Depth Sensitivity Depends on NA and Focal Plane Depth when the Functionally Activated Column is Small and Independent of NA and Focal Plane Depth when the Activated Column is Large

Here, we examine how the NA and focal plane depth affect the depth sensitivity for absorption (Figs. 8 and 9) and fluorescence imaging (Figs. 10 and 11). Figure 8 depicts the depth sensitivity for functionally activated columns of four different cross-sectional areas [Fig. 8a: 10 × 10, Fig. 8b: 50 × 50, Fig. 8c: 100 × 100, and Fig. 8d: 200 × 200 μm^{2}]. In Figs.
8a, 8b, 8c, 8d, NA is fixed at 0.2 and focal plane depth varies among 100 (circle), 300 (star), and 500 μm (triangle). Here, we plot the percent contribution of each layer to the total signal at the detector center using Eqs. 8, 9. For example, for absorption imaging with NA of 0.2 and focal plane depth of 100 μm, the top 100 μm (0–100 μm) of the tissue contributes to ~85% of the total signal (contributions from 0–1000 μm). In Fig. 8a, where the activated column represents a single array of point absorbers [see Fig. 2a] (in our simulation, the horizontal size of one voxel is 10 μm), the depth sensitivity varies drastically with the focal plane depth. In particular, the increase of the focal plane depth from 100 to 300 μm results in a significant increase in the percent contribution from the deeper tissues; the increase of the focal plane depth from 300 to 500 μm, however, results in a smaller increase of the percent contribution from the deeper tissues. For example, the contributions of the tissue within the top 100 μm are 84, 46, and 37%, while those from 200 to 300 μm are 2, 10, and 16%, corresponding to focal plane depths of 100, 300, and 500 μm, respectively. This trend of increasing percent contribution from deeper layers with the focal plane depth reaches a plateau for focal plane depth of >500 μm (not plotted here).

As the cross-sectional area of the activated column increases, the depth sensitivity still varies, but less drastically, with the focal plane depth [Figs. 8b and 8c]. When the area of the activated column reaches 200 × 200 μm^{2} [Fig. 8d], the depth sensitivity barely depends on the focal plane depth. Note that the horizontal size 200 μm is larger than the spatial resolutions for the given NA and focal plane depths. As can be seen from all panels, >97% of the signal comes from the top 500 μm of tissue; thus, only the depth sensitivity from the top 500 μm is plotted in all panels and the Figs.
8, 9, 10, 11.

Figure 9 depicts the depth sensitivity for functionally activated columns with different cross sectional areas [Fig. 9a: 10 × 10, Fig. 9b: 50 × 50, Fig. 9c: 100 × 100, and Fig. 9d: 200 × 200 μm^{2}), when the focal plane depth is fixed at 300 μm and the NA varies between 0.1 (circle), 0.2 (star), and 0.4 (triangle). The results are similar to those of Fig. 8. In Fig. 9a, where the activated column represents a single array of point absorbers, the depth sensitivity varies significantly with the NA. In particular, the increase of the NA from 0.1 to 0.2 results in a large increase of the percent contribution from the deeper tissues; the increase of the NA from 0.2 to 0.4, however, results in a smaller increase of the percent contribution from the deeper tissues. This trend of increasing percent contribution from deeper layers with the NA reaches a plateau for NA ≥ 0.4 (not plotted here). The depth sensitivity for the column with a cross-sectional area of 50 × 50 μm^{2} [Fig. 9b] is similar to that of [Fig. 9a]. When the area of the activated column increases to 100 × 100 μm^{2}, the depth sensitivity still varies, but less drastically, with the NA. When the area of the activated column reaches 200 × 200 μm^{2} [Fig. 9d], the depth sensitivity barely depends on the NA.

In summary, Figs. 8 and 9 demonstrate that the smaller the cross-sectional area of the activated column is, the more sensitive the depth sensitivity is to the NA and focal plane depth. Interestingly, the depth sensitivity does not depend on the NA or focal plane depth when the horizontal size of an activated column is larger than its spatial resolution, which can be explained as follows: δ*I*
_{z}
_{0}(0, 0, *z*
_{d}), the individual contribution to the intensity at the detector center, is proportional to the summation of the Green function *G*(*x*, *y*, *z*|0, 0, *z*
_{d}) over the cross-sectional area of the activated column. Because this Green function represents the light distribution inside the tissue from a point source at the center of the detector, its summation over a cross-sectional area with a horizontal size larger than spatial resolution is approximately the total power transmitted through the tissue, whose profile versus depth (that is, the depth sensitivity) does not change with either NA or focal plane depth.

Next, we present the results for fluorescence imaging. For fluorescence imaging with the uniform distribution of fluorescent targets, these are similar to the results shown in Figs. 8 and 9. Therefore, we only show the results for the nonuniform distribution of fluorescent targets [as depicted in Fig. 2d]. Here, we observe the same trend of the depth sensitivity change with the change of focal plane depth (Fig. 10) and NA (Fig. 11). In addition, there is a notable difference between fluorescence imaging of tissues with uniform and nonuniform staining profile. Figure 10d depicts a higher percent contribution (53%) to the signal from deeper tissues (200–500 μm) than that shown in Fig. 8d (26%), a direct consequence of the different spatial distributions of the fluorescent targets in these two cases, as shown in Figs. 2c and 2d. This illustrates the strong effect of the dye-concentration profile on the depth sensitivity. In many functional fluorescence imaging experiments, the dyes are externally administered. Hence, different dye concentration profiles may be obtained depending on the specific experimental protocols.^{20, 21, 22}

## 3.3.

### Larger NA and Focal Plane Depth Suppress the Contribution from Surface Vessels in the Absorption Imaging

As we have seen from Figs. 8d and 10d, the depth sensitivity does not depend on the NA and focal plane depth when the activated column is large. This is only true, however, when the activated columns are devoid of surface vessels above them. In this case, we can assume a uniform distribution of absorption coefficient in the tissue. However, there are areas in the brain tissue that are covered with surface vessels, which can significantly affect the absorption coefficient of that area. Here, we consider only the situation in which the horizontal size of the activated column is larger than the spatial resolution: in particular, in which it is 400 μm. The vessel of 50 × 50 × 50 μm^{3} is embedded within the top 400 μm of the tissue, and a blood tissue volume ratio of ~1% is used. Hence, the extra point absorbers due to the presence of the vessel is

## Eq. 19

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta \mu _{\rm a} (r) = \sum\limits_{\scriptstyle - 25 \le x_i \le 25\,\,{\rm \mu m} \hfill \atop {\scriptstyle - 25 \le y_j \le 25\,\,{\rm \mu m} \hfill \atop \scriptstyle 0 \le z_k \le 50\,\,{\rm \mu m} \hfill}} {100\Delta } \mu _{{\rm a}0} \delta (x - x_i, y - y_j, z - z_k). \end{equation}\end{document} $$\delta {\mu}_{\mathrm{a}}\left(r\right)=\sum _{\genfrac{}{}{0pt}{}{{\scriptstyle -25\le {x}_{i}\le 25\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\mu \mathrm{m}}}{\genfrac{}{}{0pt}{}{{\scriptstyle -25\le {y}_{j}\le 25\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\mu \mathrm{m}}}{{\scriptstyle 0\le {z}_{k}\le 50\phantom{\rule{0.16em}{0ex}}\phantom{\rule{0.16em}{0ex}}\mu \mathrm{m}}}}}100\Delta {\mu}_{\mathrm{a}0}\delta (x-{x}_{i},y-{y}_{j},z-{z}_{k}).$$Figure 12a depicts the percent contribution from different depths of tissue and the surface vessel in the absorption imaging for NA of 0.2 and focal plane depths of 100 (circle), 300 (star), and 500 μm (triangle). As we focus deeper into the tissue, the contribution of the surface vessels decreases. For example, the contribution from the surface vessel is suppressed from ~36 to 14% as the focal plane depth changes from 100 to 500 μm. A similar trend is observed for focal plane depth = 300 μm and NA of 0.1 (circle), 0.2 (star), and 0.4 (triangle), as shown in Fig. 12b. Note that when the surface vessels are either very large (≥100 μm) or small (~10 μm), the contribution of the surface vessels to the signal will be either too large or too small; thus, the effect of changing NA and focal plane depth will be negligible.

Here, we consider absorption imaging only because in many fluorescence imaging experiments the wavelength is chosen so that the absorption of the hemoglobin, and thus, of the surface vessels, is negligible.

## 4.

## Discussion

## 4.1.

### Method validation

The Monte Carlo code tMCimg itself has been validated against an accepted analytic solution for a semi-infinite medium, and cross-validated for a slab geometry with an absorbing inclusion.^{19}

We evoke the first Born approximation in our calculation. In order for this approximation to be valid, the spatial perturbation to the absorption coefficient, δμ_{a}(r), needs to be much smaller than the background absorption coefficient μ_{a0}.^{16} As discussed above, δμ_{a}(*r*) is often one order of magnitude smaller than μ_{a}(*r*); thus, the first Born approximation is valid and δ*I* in Eq. 1 can truly represent the difference image. When the perturbation is large, however, a full Monte Carlo simulation needs to be employed. In order to reduce the number of Monte Carlo simulations for calculating the spatial resolution, we take advantage of the 4*f* configuration so as to simplify Eqs. 5, 10 into Eqs. 6, 14. Here, we examine what will happen if the imaging system does not have a 4*f* configuration. In general, the two lenses L_{1} and L_{2} can have different focal lengths *f*
_{1} and *f*
_{2}, respectively and are separated by a distance *d*. Let us examine two situations: (*i*) *d = f*
_{1} *+ f*
_{2} and 2) *d ≠ f*
_{1} *+ f*
_{2}. When *d = f*
_{1} *+ f*
_{2}, the imaging system can be viewed as a general 4*f* configuration with magnification *M = f*
_{2}/*f*
_{1}. The spatial profiles and spatial resolutions shown in Figs. 6 and 7 are valid as long as we understand that they apply directly to the brain tissue. For example, if *M* = 2, then a spatial resolution of 200 μm means that we can resolve two cortical columns separated by 200 μm, though on the detector plane they will actually be separated by 400 μm (200 μm × 2). When *d ≠ f*
_{1} *+ f*
_{2}, the spatial resolution provides good approximation. For example, the spatial resolution calculated by Polimeni
^{9} using a tissue model that is similar to a particular configuration used in our simulation, though with a notable difference: in their model, *f*
_{1} *= f*
_{2}, but *d* ≪ *f*
_{1} *+ f*
_{2}, which represents a case very different from the 4*f* configuration. However, their result of 240 μm is similar to our result of 200 μm. In addition, Orbach and Cohen^{10} reported an experimental result of 280 μm using a condition similar to that of Polimeni.^{9} Hence, our simulation improves the efficiency while maintaining results that agree with other calculations and experimental measurements. Note that we do not need to assume a *4f* configuration in deriving the equations for calculating depth sensitivity; thus, the depth sensitivity for other optical configurations will be the same as shown in Figs.
8, 9, 10, 11, 12.

## 4.2.

### Computational Efficiency

We have simulated 21 different optical configurations—three values for NA (0.1–0.4) and seven for focal plane depth (20–600 μm)—covering the typical range encountered in camera-based optical brain-imaging experiments. We have considered the optical disturbance within the top 1 mm of the brain tissue and have chosen a voxel size of 10 μm. Therefore, we must calculate the contributions from 100 locations in order to obtain the spatial resolution.

As discussed above, in order to obtain the difference images, we calculate the Green function originating from the center of the detector toward the tissue [TeX:] $G(x,y,z,\hat \Omega |0,0,z_{\rm d})$ $G(x,y,z,\widehat{\Omega}|0,0,{z}_{\mathrm{d}})$ instead of that originating from the point absorbers (or fluorescent targets) inside the tissue toward the detector [TeX:] $G(x_{\rm d}, y_{\rm d}, z_{\rm d} |0,0,z,\hat \Omega)$ $G({x}_{\mathrm{d}},{y}_{\mathrm{d}},{z}_{\mathrm{d}}|0,0,z,\widehat{\Omega})$ . Here, we will compare the number of simulations needed for both approaches. In the previous case, we need one simulation for each combination of NA and focal plane depth; however, the contributions from different locations inside the tissue are automatically taken care of by [TeX:] $G(x,y,z,\hat \Omega |0,0,z_{\rm d})$ $G(x,y,z,\widehat{\Omega}|0,0,{z}_{\mathrm{d}})$ . Hence, we need a total of 21 simulations for estimating the spatial resolution and depth sensitivity of both the absorption and fluorescence imaging. In the latter case, we need 100 simulations. By using the former approach, we improve the efficiency by five times.

## 4.3.

### Strategies to Increase the Penetration Depth of Absorption and Fluorescence Imaging

As discussed in Sec. 3, changing NA and focal plane depth does not affect depth sensitivity when the functionally activated cortical areas are large and devoid of large surface vessels. However, as can be seen in Fig. 12, when imaging cortical areas covered with large surface vessels, focusing deeper into the tissue and/or enlarging the NA may suppress the contributions from these vessels and effectively enhance the penetration depth of the absorption imaging. In fact, this practice has been employed in functional imaging to reduce the effect of the surface vessels, as demonstrated in Ref. 15.

In many camera-based optical brain imaging experiments, large cortical areas are often excited during the peak response. For example, electrically stimulating the forepaw of a rat can lead to functional activation of a cortical area well beyond 1 mm^{2}. Hence, changing the NA and focal plane depth does not affect the penetration depth (Fig. 8). On the other hand, with the advent of new technologies, such as optogenetics in which optical methods are used to control genetically targeted neurons within intact neural circuits, it is likely that selective activation of a small cross-sectional area ≤ 100 × 100 μm^{2} may be soon achieved.^{26} As shown in Figs.
8, 9, 10, 11, focusing deeper into the tissue and/or enlarging the NA would enhance the penetration depth of both absorption and fluorescence imaging of such small cortical areas.

In fluorescence imaging, where dyes are applied externally, it may be possible to enhance the contributions from specific cortical layers by regulating the density distribution of the dyes across the cortical depth. For example, in voltage sensitive dye imaging, Ferezou
^{20} reported a density distribution shown in Fig. 2d, where deeper layers (200–500 μm) had more dye than the shallow layers (0–200 μm). Therefore, the image would consist of more contribution from deeper layers [comparing the two profiles of the depth sensitivity, in Figs. 8d and 10d].^{20}

## 4.4.

### Depth Sensitivity is Important in Data Interpretation, Cross Modality and Species Comparison of Functional Imaging Results

It is common to compare functional brain-imaging data from different imaging modalities (such as IOI and VSDI) or across species (such as mouse and rat). In general, the mammalian neo-cortex consists of six layers, each with a unique thickness depending on the species and cortical areas.^{27} To compare fairly, data collected from different modalities or species should come from the same cortical layer(s). Hence, it is critical to know quantitatively the contribution from each layer with a specific imaging modality (e.g., IOI) on a particular animal model (e.g., mouse). Here, we discuss two examples and determine for each whether it is appropriate to compare data obtained from two different modalities or species.

We first compare the data from the primary somato-sensory cortex (SI) of a rat obtained from two different modalities: IOI and VSDI. For simplicity, we only consider the case when the activated column is larger than or equal to the spatial resolution and the staining profile for VSDI is nonuniform, as shown in Fig. 2d. As discussed earlier, 97% of the contribution to the signal comes from the top 500 μm. Therefore, we need to consider only layer I (0–150 μm) and II/III (150–650 μm) of the rat SI.^{27} By using Eqs. 8, 9 for IOI and Eqs. 17, 18 for VSDI, we have obtained the percent contribution from layer I and II/III to the total signal to be 58 and 41% when using IOI, and 26 and 73% when using VSDI. Hence, although the signal is dominated by layer I in IOI, it is dominated by layer II/III in VSDI.

The second example examines the IOI results obtained from the same area of two different species: the SI of a mouse and a rat. For simplicity, we only consider the case when the activated column is larger than the spatial resolution. The quantitative analysis of the depth sensitivity shows that 38 and 56% of the signal are from layer I (0–100 μm) and II/III (100–400 μm) of the mouse SI, respectively.^{27} Once again, although we see that the signal is dominated by layer I in the rat, it is dominated by layer II/III in the mouse.

Both examples reveal that caution has to be taken in comparing data across modalities and species. In general, we need to estimate the contribution from each cortical layer of an individual case in order to compare meaningfully data across different imaging modality and species.

## 4.5.

### Knowledge of Spatial Resolution and Depth Sensitivity Helps Guide Optical Imaging Optimization

When choosing a camera-based optical imaging system, we need to consider important parameters such as NA, focal plane depth, and pixel size of the detector. Our results (Figs. 7, 8, 9, 10, 11, 12) provide both qualitative and quantitative guidance as to how to optimize these parameters: the increase of NA and focal plane depth results in a lower spatial resolution, more suppression of the contribution from the surface vessels, and more contribution from deeper tissues (when the activated column is sufficiently small). In general, larger NA also leads to a stronger signal. Here, we consider some typical examples and give quantitative estimation for each case, as follows:

*Imaging a cortical region covered by surface vessels.*We can increase NA or focal plane depth, or both, when we want to minimize the contribution of surface vessels and when the spatial resolution is not of primary concern. For example, in imaging the cortical region corresponding to the rat forepaw stimulus, we can use NA = 0.4 and focal plane depth = 300 μm to effectively suppress the contribution of the surface vessel while maintaining a spatial resolution of ~200 μm.*Imaging a cortical region devoid of surface vessels.*Here, if a fine spatial resolution is required to differentiate two neighboring cortical columns such as two barrels, we may want to use NA = 0.4 and focal plane depth = 100 μm to maximize the signal strength and keep the spatial resolution to ~ 50 μm, sufficient to resolve two columns. After we know the spatial resolution of the system, we can specify the largest pixel size according to the Nyquist sampling theorem. For example, if the spatial resolution is 200 μm and the magnification of the imaging system is 1, the largest pixel size is 100 μm. This is important in selecting the right type of camera and in optimizing its speed and resolution.

## 5.

## Summary

We have estimated the spatial resolution and the depth sensitivity of the camera-based 2-D optical imaging methods for a wide range of optical configurations (NA: 0.1–0.4; focal plane depth: 20–600 μm). We have found the following: (*i*) the spatial resolution worsens with increases in NA and focal plane depth; (*ii*) for NA ≤ 0.2 or focal plane depth ≤ 300 μm, the spatial resolution is <200 μm, sufficient to resolve individual functional columns; (*iii*) for functional imaging applications with large activated columns devoid of surface vessels, increasing NA and focal plane depth does not increase the contribution from deeper tissues, but rather reduces the spatial resolution; (*iv*) increasing NA and focal plane depth may improve the contribution from deeper tissues when small columns are activated or when activated columns are covered with large surface vessels; and (*v*) our results offer guidance in data interpretation and cross-species/modality comparison, as well as in the optimization of optical imaging configurations.

## Acknowledgments

We thank Ivan Teng and Dr. Qianqian Fang for their help on the simulation program. We gratefully acknowledge support from the NINDS (Grant Nos. NS-051188 and NS-057198 to A. D. and Grant No. NS-057476 to D. A. B.) and NIBIB (Grant No. EB-009118 to A. MD.). We also acknowledge the support from John Carroll University through the Summer Research Fellowship to P. T.

## References

*In-vivo*Optical imaging of cortical architecture and dynamics,” Modern Techniques in Neuroscience Research, 893 –969 Springer, New York (1999). Google Scholar

*in vivo*: techniques and applications from animal to man,” J. Biomed. Opt., 12 (5), 051402 (2007). https://doi.org/10.1117/1.2789693 Google Scholar