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 Koester1 for imaging the cornea; and the theta microscope, developed by Stelzer and Lindek2 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 melanin7, 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.
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 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 of RAM. The queue system is a Lava/Sun Grid Engine, running the ROCKS OS. Average simulation time was and .
The computational area is the part of the model that holds both the structure and index of refraction values that are used in the simulation. The computational area is built in steps, starting with a homogeneous body of water with . Then, additional components are inserted into the computational area, as shown in Fig. 2 . A single bead with at a depth of , as seen in Fig. 2a, is used as a test target. The radius of the bead used throughout all the experiments was . The Airy disk radius was calculated asis the free-space wavelength, ; is the index of refraction of water, 1.33; is the diameter of the capture side of the pupil, ; and is the focal length, . The bead was chosen because a bead smaller than the Airy disk would behave as a point scatterer, and it did not produce sufficient signal at the detector. The use of a smaller particle would have resulted in signals too small to evaluate in comparison to the background arising from light scattered by out-of-focus objects. A larger particle would have provided stronger signals, but as shown later, it is not always easy to predict where to focus the microscope to achieve the strongest signal. A flat DE junction is created by filling the computational area below with a medium that has in Fig. 2b. For a more realistic model of the DE junction a sinusoidal junction [Fig. 2c] is built at a mean depth of with a amplitude of and a period of , and . Next, cells and some subcellular components were added, as shown in Fig. 2d, which has a uniform dermis with , an epidermis built with cytoplasm with , a cell nucleus with , and the intercellular fluid with (Refs. 9, 10, 11). Figure 2e includes mitochondria, which have (Ref. 12). Figure 2f includes the melanin with (Ref. 13), and is the most detailed skin model we use. Both the mitochondria and melanin are placed pseudorandomly for each row of cells.
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 radius and vertical minor axis with a radius of (Ref. 14). Cells close to the dermis are circles14 of radius . The nucleus of each cell is 10% of the total cell volume.15 The clumps of pigment protein melanin have radii of 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 and minor axis radii of 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 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 is perpendicular to the plane, the component points down along the ordinate, and is toward the right along the abscissa. Then, this computational area was used in the FDTD script to com-pute the , , and 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 , 21 focal locations were chosen. Each focal location was 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 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, from the last bead, and the simulation is run again. The shifting process is illustrated in Fig. 4 .
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 in a uniform medium of water with an 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.
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 in the direction and a width of , or 33% of the width of the computational area. The amplitude of the source is2 is modulated by the ramp function18 2, is the location in the lateral direction in the pupil plane; is the vertical location of the pupil; is the current time; the angular frequency is 2.69 peta-radians per second, which corresponds to a vacuum wavelength of ; is , which represents the farthest location laterally in the computational area; is the wave number; and is the focal length of the Gaussian wave. In Eq. 3, is and equals the time when the error function (erf) is at the inflection point, and the rise time has a value of . Using the ramp function [Eq. 3], the simulation reaches steady state sooner than it would if it was required to wait for the starting transient to decay. The Gaussian is defined only in the left half of the pupil. The location of the source and receiver is from the top of the computational area.
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 mode wave propagation using the discretized forms of the following lossless curl equations:is the permittivity of the medium, and is the permeability of the medium. In Eqs. 7, 8, 9, is the time step index, ( and ) are the coordinates in the computational area, is the temporal step size, and is the spatial step size. For the FDTD to remain stable, two conditions must be met. First the Courant condition must be satisfied such that is the temporal step size; is the spatial step size; the comes from the number of dimensions used in the simulation; is the speed of light in free space; and is the minimum index of refraction, 1.33, which gives the smallest time step.20, 21 The second condition for stability in the FDTD is that the spatial step size must be less than the wavelength , generally by a factor of 10, leading to is 1.7, the highest value in the computational area; and is the free-space wavelength, (Ref. 21). Thus, Eq. 10 becomes of (Refs. 22, 23). Depending on the source location, this simulation used either a perfectly matched layer (PML) or a Mür first-order absorbing boundary condition22, 24, 25 (ABC). When the excitation source was within of the absorbing boundary (AB), the Mür first-order ABC was used instead of the PML, as a PML is not designed to handle24 sources close to the AB.
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 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 of the microscope. For example the theta line-scanning microscope has a of , so a pinhole has an image at the sample compared to the ideal Airy disk diameter of . 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,is is the focal length, as in Eq. 2; is the location in the pinhole; and is the location in the array in the pupil. Integrating Eq. 13 using Simpson’s rule yields is is the location in the array in the pupil, as in Eq. 14. To generate the diffraction pattern at the image plane, Eq. 15 is applied to each location . The squared absolute value of the result from Eq. 15 is the irradiance. Numerically integrating over all points subtended by the pinhole results in the detected power for a specific pinhole size.
Imaging Performance Measure
In the simulation presented here, the signal is defined as the detected power when there is a bead at the focal point, and the clutter as that when there is no bead at the focal point. As 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 of 1.33. This noise determines the lower limit of the simulation’s detection range.
The signal-to-clutter ratio is defined asis received signal power for a specific pinhole diameter, and is the received clutter power signal for a specific pinhole diameter. Note that includes the contribution from clutter , 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.
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.
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 . Then, the bead was moved in steps with maximum distances of and from , along the ordinate axis, and a maximum distance of 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 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
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 , 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 , as shown in Fig. 7 . In all of these examples, with one exception, for mitochondria, the was higher than . Figure 5 plots the and separately for the mitochondria test case. This figure shows greater than and confirms that the was equal to or lower than the for pinhole widths near . 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 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 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 pinhole (shown in Fig. 9 ) we see there are 10 locations where the 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 . 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 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 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.
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 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.
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.