## 1.

## Introduction

Finite-difference time domain (FDTD) simulations are useful tools in exploring how light interacts with complex scattering media. In this paper, we used an FDTD simulation to explore how light from a confocal reflectance theta line-scanning microscope travels from the source to the focus, and then to the receiver. The design of the confocal reflectance theta line-scanning microscope is based on two existing microscope designs: the line-scanning microscope with a divided objective lens, developed by Koester^{1} for imaging the cornea; and the theta microscope, developed by Stelzer and Lindek^{2} and Webb.^{3} Theta line-scanning microscopes differ from conventional confocal microscopes by being confocal in one direction and nonconfocal in the other. The theta line-scanning microscope may lead to a smaller, more clinically usable confocal reflectance microscope for dermatological and other applications. One way this device decreases the size and cost is by changing the usual point source and point detector scheme to a line source with a 1-D array detector. Thus, this configuration requires only a 1-D scan, thereby producing confocality in one direction, but not in the other. Figure 1
is a schematic of the light path for a recently developed theta line-scanning microscope.^{4, 5} This figure shows the bistatic nature of the microscope, in which the source illumination and detected wave traverse different paths through the sample.

In comparison with the conventional confocal microscope, the theta line-scanning microscope has increased variations in the signal, especially below the dermal-epidermal (DE) junction. These abnormalities in the received signal are displayed as darkened areas and in extreme cases, the received signal may be below the noise, thus undetectable. The darkened areas are localized to specific locations within the sample, but are larger than isolated pixels, so their origin is unlikely to be speckle. We hypothesize that these variations are caused by differences in the propagation paths of incident and received scattered light through different cells and cell components, and differing angles of incidence along the somewhat sinusoidal DE junction.

The computer model developed to test this hypothesis is a FDTD computational model. Past computer models developed to examine confocal microscopes were generally limited to Monte Carlo models due to computational constraints.^{6} On one hand, Monte Carlo models rely on averaging over a number of scattering events, which is advantageous in the sense that they enable results to be obtained with only low-order statistical information about the medium of propagation. On the other hand, they fail to model the refraction at cell boundaries, index gradients, and the multiple sizes and shapes of tissue components. They work well when most of the light detected has scattered many times. FDTD is a more realistic model, and enables us to model the signal from a scattering particle seen through a few cells with any desired degree of fidelity. However, to exercise this ability, the model requires more information, such as specific locations of different materials, and this information is often only known statistically. We believe that there is a role for both types of models, and that FDTD provides valuable insights into the behavior of confocal microscopes.

The computational model, in contrast to physical experiments, facilitates introduction of skin components individually, or in groups, to test the effect of each component on the received signal. Specifically, the model revealed that, although different cell components and the DE junction all contribute to the received signal reduction, the presence of the pigment protein melanin^{7, 8} causes the greatest reduction. In addition, we explored the scattering of light from a single sphere using the model, with results different than those of a conventional, common-aperture confocal reflectance microscope.

## 2.

## Methods

The simulation was constructed in four major computational blocks. The first block builds the computational area, containing the physical description of the media being modeled. Next is the actual FDTD kernel, which is used to propagate the electromagnetic wave in the computational area. Following the FDTD kernel is the Fresnel-Kirchhoff integral, used in calculating the irradiance at the pinhole plane. The last computational block deals with the calculation the power received at the detector.

The simulations were completed in MATLAB version 7.0. Most of the simulations were performed on a Fedora Core 4 $64\text{-}\text{bit}$ PC running the Linux 2.6 kernel. The simulation of the scattering of light from a sphere required several hundred simulations. These were performed on a Linux distributed memory cluster named “Opportunity” at Northeastern University. The cluster has 64 nodes, each with two CPUs and $4\phantom{\rule{0.3em}{0ex}}\mathrm{Gbytes}$ of RAM. The queue system is a Lava/Sun Grid Engine, running the ROCKS OS. Average simulation time was $2\phantom{\rule{0.3em}{0ex}}\mathrm{h}$ and $40\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ .

## 2.1.

### Computational Area

The computational area is the part of the model that holds both the structure and index of refraction $\left(n\right)$ values that are used in the simulation. The computational area is built in steps, starting with a homogeneous body of water with $n=1.33$ . Then, additional components are inserted into the computational area, as shown in Fig. 2 . A single bead with $n=1.47$ at a depth of $45\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , as seen in Fig. 2a, is used as a test target. The radius of the bead used throughout all the experiments was $0.8\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ . The Airy disk radius was calculated as

## Eq. 1

$${R}_{\text{Airy}}\approx \left(\frac{1}{2}\right)\left(1.22\right)\frac{\lambda}{nD(0.5\u2215f)}\approx 0.577\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m},$$The sizes of the cells and cell components in the skin model are chosen to match their typical physical dimensions in human skin cells. Cells close to the surface of the epidermis are elliptical, with a horizontal major axis of approximately
$20\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
radius and vertical minor axis with a radius of
$5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
(Ref. 14). Cells close to the dermis are circles^{14} of radius
$10\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
. The nucleus of each cell is 10% of the total cell volume.^{15} The clumps of pigment protein melanin have radii of
$0.4\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
and occupy in total 8.5% of the cell volume.^{15} The melanin clumps are the smallest element in the simulation. The mitochondria are also approximated by ellipsoids, with major axis radii of
$2\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
and minor axis radii of
$0.75\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
in the transverse direction, and occupy 10% of the total cell volume.^{15, 16, 17} In our 2-D model we used the given volume ratios for an ellipse to find the volume ratios for a cylinder. From the cylinder volume ratios we calculated the 2-D area ratios.

In the fully modeled skin in Fig. 3 , we assume the medium above the epidermis has an $n$ of 1.33. In practice, an index-matching gel is used with approximately the same index. Additional structures were not included in the dermis at this time because we are currently interested in imaging near the DE junction. In the computational area in Fig. 3, the ${E}_{z}$ is perpendicular to the plane, the ${H}_{y}$ component points down along the ordinate, and ${H}_{x}$ is toward the right along the abscissa. Then, this computational area was used in the FDTD script to com-pute the ${E}_{z}$ , ${H}_{x}$ , and ${H}_{y}$ fields.

To mimic the scanning process of the theta line-scanning microscope a variety of different focal locations were simulated. The computational area of Fig. 2f was used for this experiment. Working at a depth of $45\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , 21 focal locations were chosen. Each focal location was $4\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ from its neighbor. These 21 focal locations approximate the scanning path of the theta line-scanning microscope. To simulate each focal location, the computational area is shifted to the left by $4\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ and the gap on the right side is filled with more cells, such that there are no discontinuities. With the gap filled, a new bead is placed at the focal location that is to be sampled, $4\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ from the last bead, and the simulation is run again. The shifting process is illustrated in Fig. 4 .

## 2.2.

### Experiments Performed

Two computational experiments were conducted with this model. The first uses a spherical bead with an index of refraction of 1.48 and a radius of $0.8\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ in a uniform medium of water with an $n$ of 1.33. With the focus location held fixed, the bead was moved to various locations in the medium to simulate the scanning effects of the microscope. The second experiment involved shifting a high-fidelity skin model 21 times to examine variations in the signal and clutter intensities.

## 2.3.

### FDTD

The upper edge of the computational area is considered the pupil of the microscope. The pupil is split, with the source on the left and the detector optics on the right. The model uses the FDTD algorithm to simulate light from the source on the left half of the pupil to the target bead and back to the receiver pupil on the right half of the pupil. The light source in the theta line-scanning microscope is modeled as a converging spherical Gaussian wave, traveling from the left half of the upper boundary of the computational area into the skin, as seen in Fig. 3. The excitation source has a peak amplitude of $1\phantom{\rule{0.3em}{0ex}}\mathrm{V}\u2215\mathrm{m}$ in the $z$ direction and a width of $34.31\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , or 33% of the width of the computational area. The amplitude of the source is

## Eq. 2

$${E}_{z}(x,{y}_{0})=\mathfrak{R}\left({f}_{\mathrm{erf}}\left(t\right)\mathrm{exp}\left(i\omega t\right)\mathrm{exp}{[-\frac{x-({X}_{\mathrm{max}}\u22154)}{(2{X}_{\mathrm{max}}\u22156)}]}^{2}\mathrm{exp}\{ik{[{f}^{2}+{(x-\frac{{X}_{\mathrm{max}}}{2})}^{2}]}^{1\u22152}-f\}\right).$$^{18}

## Eq. 3

$${f}_{\mathrm{erf}}\left(t\right)=\frac{1}{2}[1+\mathrm{erf}\left(\frac{t-{t}_{0}}{T\sqrt{2}}\right)].$$The FDTD algorithm was developed by Yee as a numerical solution to the time-dependent Maxwell’s equations.^{19} Yee’s FDTD can be used to solve most generic wave propagation problems in three dimensions. This work does not use the full 3-D FDTD solution to Maxwell’s equations, but rather, the simpler 2-D FDTD solution. The simulation models the 2-D
${\mathrm{TM}}_{z}$
mode wave propagation using the discretized forms of the following lossless curl equations:

## Eq. 4

$$(\frac{\partial {H}_{y}}{\partial x}-\frac{\partial {H}_{x}}{\partial y})=\u03f5\frac{\partial {E}_{z}}{\partial t},$$## Eq. 5

$$\left(\frac{\partial {E}_{z}}{\partial y}\right)=-\mu \frac{\partial {H}_{x}}{\partial t},$$## Eq. 7

$${E}_{z}^{m+1}(i,j)={E}_{z}^{m}(i,j)+\frac{\Delta t}{\u03f5(i,j)\Delta}\times [{H}_{y}^{m+1\u22152}(i+\frac{1}{2},j)-{H}_{y}^{m+1\u22152}(i-\frac{1}{2},j)+{H}_{x}^{m+1\u22152}(i,j-\frac{1}{2})-{H}_{x}^{m+1\u22152}(i,j+\frac{1}{2})],$$## Eq. 8

$${H}_{x}^{m+1\u22152}(i,j+\frac{1}{2})={H}_{x}^{m-1\u22152}(i,j+\frac{1}{2})-\frac{\Delta t}{\mu \Delta}[{E}_{z}^{m}(i,j+1)-{E}_{z}^{m}(i,j)],$$## Eq. 9

$${H}_{y}^{m+1\u22152}(i+\frac{1}{2},j)={H}_{y}^{m-1\u22152}(i+\frac{1}{2},j)+\frac{\Delta t}{\mu \Delta}[{E}_{z}^{m}(i+1,j)-{E}_{z}^{m}(i,j)],$$^{20, 21}The second condition for stability in the FDTD is that the spatial step size must be less than the wavelength $\lambda $ , generally by a factor of 10, leading to

## Eq. 11

$$\Delta =\frac{\lambda}{{n}_{\mathrm{max}}10}=41.17\phantom{\rule{0.3em}{0ex}}\mathrm{nm},$$## Eq. 12

$$\Delta t\u2a7d\frac{\lambda}{10\sqrt{2}{c}_{0}({n}_{\mathrm{max}}\u2215{n}_{\mathrm{min}})}=0.129\phantom{\rule{0.2em}{0ex}}\mathrm{fs}.$$^{22, 24, 25}(ABC). When the excitation source was within $411\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ of the absorbing boundary (AB), the Mür first-order ABC was used instead of the PML, as a PML is not designed to handle

^{24}sources close to the AB.

## 2.4.

### Fresnel-Kirchhoff Integral

When the received backscattered light travels through a lens to the image plane, the diffraction pattern must be calculated using the Fresnel-Kirchhoff integral. Because the Fresnel-Kirchhoff integral uses a complex field, the time domain (TD) signal from the FDTD simulation was converted to a frequency domain (FD) signal. During the simulation’s last time period near steady state, the ${E}_{z}$ electromagnetic field in the pupil plane was converted using the Fourier transform to obtain the FD signal data. The FD signal data were then propagated backward through a uniform, isotropic medium to the image plane. Because the pinhole is conjugate to the transmitter focal point, propagating the wave backward is mathematically identical to propagating the wave through a lens and focusing it onto a pinhole, as shown in Fig. 3. The diameter of the pinhole used in the simulation was determined by dividing the diameter of the physical pinhole by the magnification $M$ of the microscope. For example the theta line-scanning microscope has a $M$ of $6\times $ , so a $50\text{-}\mu \mathrm{m}$ pinhole has an $8.33\text{-}\mu \mathrm{m}$ image at the sample compared to the ideal Airy disk diameter of $1.154\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ . The irradiance at the pinhole was computed and plotted to show the effects of propagation through the medium.

The Fresnel-Kirchhoff integral in two dimensions,

## Eq. 13

$$E\left(x\right)=\frac{2ki}{4\pi}{\int}_{\text{pupil}}{E}_{z}\left(x\right){e}^{ikr}\phantom{\rule{0.2em}{0ex}}\mathrm{d}x,$$## Eq. 15

$$E\left(x\right)\simeq \frac{2ki}{4\pi}\sum _{q=1}^{m}{E}_{z}\left[\frac{\Delta}{3}({s}_{q-1}+{s}_{q}+{s}_{q+1})\right],$$## Eq. 16

$${s}_{q}=\mathrm{exp}\left\{ik{[{(x\left(q\right)-{x}_{p})}^{2}+{f}^{2}]}^{1\u22152}\right\}.$$## 2.5.

### Imaging Performance Measure

In the simulation presented here, the signal
${S}_{p}$
is defined as the detected power when there is a bead at the focal point, and the clutter
${C}_{P}$
as that when there is no bead at the focal point. As
${C}_{P}$
arises from scatterers within the medium, the clutter is not the same as noise. Noise is defined in experimental systems as the variation in power measured when the source is turned off, and determines the detection statistics. There are three causes of noise in a conventional microscope: electronic noise, quantum noise, and background noise. The limiting factor in confocal microscopes is either the electronic noise or quantum noise.^{26}

Computational noise is defined in this paper as the power that is detected with a uniform computational area having an $n$ of 1.33. This noise determines the lower limit of the simulation’s detection range.

The signal-to-clutter ratio $\left({\mathrm{SCR}}_{P}\right)$ is defined as

where ${S}_{P}$ is received signal power for a specific pinhole diameter, and ${C}_{P}$ is the received clutter power signal for a specific pinhole diameter. Note that ${S}_{P}$ includes the contribution from clutter ${C}_{P}$ , but through coherent addition, so the numerator in Eq. 17 can be negative if the bead signal is small and out of phase with the clutter. Equation 17 calculates just the ratio between the received signal power excluding clutter and the clutter. Figure 5 shows the signal and clutter power at the pinhole as functions of the pinhole size for the configuration in Fig. 2e. In general, both the signal and the clutter increase as the pinhole size increases. In this particular example, the signal is always sufficient to produce a result greater than the clutter.## 3.

## Results

Results from the two computational experiments are discussed in detail here. The first set of experiments examined the received signal from a single bead in various media. The second set of experiments involved scanning the bead through the focus of the beam.

## 3.1.

### Scattering from a Bead

For the first experiment, the computational area of the simulation is similar to Fig. 2a. The bead started at a base depth of
$45\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
. Then, the bead was moved in
$0.206\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
steps with maximum distances of
$+2.47$
and
$-1.64\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
from
$45\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
, along the ordinate axis, and a maximum distance of
$\pm 1.23\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
from the center of the computational area along the abscissa. Much larger steps would have produced insufficient detail and smaller ones would require excess amounts of computational time. The power magnitude
${S}_{P}$
was plotted in two dimensions as a function of the bead’s center location, as seen in Fig. 6
. Figure 6 illustrates a distinctive feature of the theta configuration: two asymmetric peaks. In a common-aperture confocal microscope, if the wavefront is not distorted, the greatest signal is expected when the wavefront of the incident wave matches the surface of the bead, i.e., the microscope is focused at the center of the bead. The reflected waves appear to originate from a point at the center of the bead and, therefore, pass through the pinhole conjugate to the laser source. If the bead is moved to any other position, the reflected wave is not matched to the incident wave, and the light is spread over a wider area in the pinhole plane, resulting in a reduced signal. Thus, there is one maximum, when the microscope is focused at the center of the bead. Figure 6 shows two received power signal peaks, one when the bead is at a shallower depth in the computational area, so that the microscope is focused at the bottom of the bead, and one when the bead is deeper, so that the focus is at the top of the bead. The received signal power peaks at the top and bottom because the surface of the bead is normal to the bisector of the transmitter and receiver. The resulting received signal power peaks come from the strong specular reflection at the surfaces of the bead. The specular reflection from the bottom of the bead is more intense due to focusing by its curved shape and additional focusing caused by the index mismatch at the convex interface at the top. Figure 6 also indicates that the detected power increases as the pinhole slit width increases, but the resolution becomes worse. The worsening of the resolution is related to the diameter of the pinhole and the width of the Gaussian beam.^{27} Generally, the pinhole is a few times larger than the transmitted Gaussian beam so that the detector receives enough light to produce sufficient signal. The unintended result of having a pinhole larger than the transmitted Gaussian beam is that the axial resolution is much worse than predicted by diffraction theory. While the transverse resolution is primarily determined by the Gaussian beam, the axial resolution is also determined by the pinhole. This is seen in Fig. 6. As the pinhole diameter is increased, the transverse resolution does not decay nearly as rapidly as the axial resolution. In particular, the optical axial sectioning decreases much faster than the lateral resolution. The simulation data for the optical axial sectioning agrees with a similar experiment performed with the theta line-scanning microscope using a mirror to determine the optical axial sectioning.^{4}

## 3.2.

### Skin Simulations

The computational area of the simulation initially contained a sole bead in a homogeneous medium, [Fig. 2a]. To understand which skin components created the lowest ${\mathrm{SCR}}_{P}$ , skin components were added to the computational area, increasing in complexity, until the full skin model was formed, as in Fig. 2f. As each additional component was added to the computational area, there was an associated drop in the ${\mathrm{SCR}}_{P}$ , as shown in Fig. 7 . In all of these examples, with one exception, for mitochondria, the ${S}_{P}$ was higher than ${C}_{P}$ . Figure 5 plots the ${C}_{P}$ and ${S}_{P}$ separately for the mitochondria test case. This figure shows ${S}_{P}$ greater than ${C}_{P}$ and confirms that the ${S}_{P}$ was equal to or lower than the ${C}_{P}$ for pinhole widths near $50\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ . These results are specific to the particular configuration of the tissue and location of the focus.

To explore the range of variations, the shifting algorithm described in Sec. 2.1 and illustrated in Fig. 4 was implemented to mimic the scanning effect of the theta line-scanning microscope. Figure 8 shows the ${\mathrm{SCR}}_{P}$ as a function of pinhole diameter for all 21 focal locations. Each curve is different from the others, as a result of the different media along the path of propagation. Note that in some cases, the ${\mathrm{SCR}}_{P}$ is negative for most of the pinhole’s diameters. The negative value indicates either that the field from the bead adds destructively to the clutter or the signal from the bead has been moved beyond the diameter of the pinhole. In such cases, the signal from that focal location cannot be detected. In other cases, the signal from the bead is detected when the diameter of the pinhole is increased. Examining the specific case of a $10\text{-}\mu \mathrm{m}$ pinhole (shown in Fig. 9 ) we see there are 10 locations where the ${\mathrm{SCR}}_{P}$ falls below zero. If the paths traveled from source to the bead to the receiver contains large-scale heterogeneous objects, then the irradiance peak is generally narrow but shifted in the pinhole plane. In these cases, a larger pinhole width enables this shifted peak to be captured, raising the ${\mathrm{SCR}}_{P}$ . If the path contains many small-scale objects, such as the melanin in this model, then the irradiance peak is generally broad, and changing the width of the pinhole in the simulation does not always cause the ${\mathrm{SCR}}_{P}$ to increase. The effect of the melanin on the irradiance peak at the detector is shown in Fig. 10 , where the irradiance peak is much broader. This results in a drop in the ${\mathrm{SCR}}_{P}$ from the melanin, as shown in Fig. 7. Because of range of scales of the heterogeneities, the optimal pinhole size fluctuates widely, as we can see in Fig. 8.

## 4.

## Conclusion

Understanding the propagation of light through skin and the effects of propagation on imaging with the theta line-scanning microscope is important to optimizing and using this technology. A 2-D FDTD simulation was used to study the propagation, as well as the imaging properties of the theta line-scanning microscope. The simulation yielded two significant results. One is that the scattering from a sole bead in a homogeneous medium results in power peaks at the top and bottom, as the bead is moved in and out of the focal point. The second was that the theta line-scanning configuration introduces extra fluctuations in signals because of the lack of a common path for incident and detected light. Extra fluctuations also came from the inclusion of melanin in the model. To some extent the effect of these fluctuations can be reduced by judicious choice of pinhole size, but they still result in regions where the signal falls below the clutter. The lower-than-expected detected signals as compared to expatiations can also be attributed to the low contrast of the bead. Increasing the index of refraction of the bead with reference to the background would have resulted in a stronger
${S}_{P}$
signal, and brought the simulation’s results closer to those images obtained by Dwyer, but the higher index would require a smaller computational spatial step size and increase computational time.^{4}

This high-fidelity 2-D model will also be useful for functional optimization of the theta line-scanning microscope, studies of adaptive optics in skin imaging, analysis of the speckle problem in reflectance confocal microscopes, and determining the signatures of different components of normal and diseased skin. Future work includes the incorporation of the measured index of refraction of real skin samples, the addition of more realistic features in the dermis, and studies of reflectance confocal microscopes at different wavelengths. Additional future work will also include simulation of a common aperture confocal microscope and comparing those results with the work presented here.

## Acknowledgments

The authors wish to thank Christopher Yocum of the University of Edinburgh for writing the Perl code for modifying the MATLAB files to run the simulations of scattering light from a sphere. We wish to thank Milind Rajadhyaksha for discussion on skin morphology. We also wish to thank Amy Moore of Massachusetts Institute of Technology and Jennifer Hogan of Tufts University for discussions concerning subcellular properties. We also thank the office of the Provost at Northeastern University for use of the “Opportunity” cluster. Professor DiMarzio’s work was supported in part by CenSSIS, the Center for Subsurface Sensing and Imaging Systems, under the Engineering Research Centers Program of the National Science Foundation (Award No. EEC-9986821). Mr. Simon was supported in part by a GAANN fellowship program entitled “Sensing, Identification, Diagnostics and Rehabilitation,” in the College of Engineering at Northeastern University.