## Abstract

The FDTD method was used to simulate focused Gaussian beam propagation through multiple inhomogeneous biological cells. To our knowledge this is the first three dimensional computational investigation of a focused beam interacting with multiple biological cells using FDTD. A parametric study was performed whereby three simulated cells were varied by organelle density, nuclear type and arrangement of internal cellular structure and the beam focus depth was varied within the cluster of cells. Of the organelle types investigated, it appears that the cell nuclei are responsible for the greatest scattering of the focused beam in the configurations studied. Additional simulations to determine the optical scattering from 27 cells were also run and compared to the three cell case. No significant degradation of two-photon lateral imaging resolution was predicted to occur within the first 40 *µ*m of imaging depth.

©2009 Optical Society of America

## 1. Introduction

Two photon fluorescence microscopy [1] has become one of the most widely used methods for fluorescence imaging due to its ability to produce images with three dimensional resolution in highly scattering samples. The two-photon technique additionally enjoys an inherent depth resolution due to the dependence of the fluorescence excitation on the square of the intensity of the incident beam [2]. Other imaging techniques such as confocal microscopy can achieve this depth discrimination in other ways, but all suffer from limits on image resolution and maximum imaging depths arising primarily from scattering. From a fluorescence excitation standpoint, scattering reduces the number of photons that reach the desired excitation volume. This results in a reduction of the intensity of the fluorescence emitted in the desired volume and can also reduce the imaging resolution of the microscope. An advantage of two-photon fluorescence over confocal techniques is that near-infrared excitation light is used, which penetrates deeper into tissue due to reduced scattering.

There has been extensive work characterizing both the point spread function (PSF) of two-photon imaging techniques and the two-photon fluorescence signals generated within the beam focal volume as a function of tissue scattering parameters [3–8]. Experimentally, it has been found that the two-photon PSF remains relatively unaffected by increased numbers of scatterers and that imaging resolution only seems to degrade in cases where the focal plane is at a depth greater than 500 *µ*m [4,5,7]. Measured fluorescence signal, however, has been shown to depend strongly on both the scattering properties of the imaging medium and depth of the focal plane [5, 7, 8].

In addition to experimental measurements, Monte Carlo models have been broadly used to characterize two-photon imaging modalities [3, 5–8]. The Monte Carlo method uses a probabilistic model to calculate the paths taken by a large number of photons in the imaging medium. The trends in fluorescence signal predicted by the Monte Carlo models generally match those found experimentally, but such models require a simplification of the problem by lumping the various optical scatterers within a tissue into bulk parameters, such as the absorbing and scattering coefficients. To our knowledge, a rigorous three-dimensional electromagnetic solution to the problem of focused optical scattering from heterogeneous biological cells has not yet been attempted.

The finite difference time domain (FDTD) method is a three-dimensional full vector computational solution to Maxwell’s curl equations [9]. While computationally costly, it is capable of modeling objects with complex physical structure and varying dielectric properties. In the scattered field formulation of FDTD [10] any incident electromagnetic field that can be formulated analytically can be easily used as an excitation source. This allows the FDTD method to be used to study the near-field scattering of a focused Gaussian beam by cells of arbitrary complex structure.

Association of particular scattering effects to specific physiological characteristics of cells has been difficult, chiefly because cells have complex structures with varying indices of refraction. In the case of far-field scattering of a plane wave by a single inhomogeneous cell, several attempts have been made to study the scattering contributions of cellular fine structure using the FDTD method [11–13]. Dunn and Richards-Kortum [11] found that organelles with diameters ranging from 0.4 *µ*m to 1 *µ*m play a significant role in the far-field scattering from a single cell, especially at angles greater than 90 degrees. Liu and Capjack [12] investigated the far-field scattering differences between large and small organelles, the effect of the cellular membrane and the role of cellular DNA content. FDTD simulations by Li et al. [14] resulted in far-field scattering signatures that were sensitive to sub-wavelength features in the scattering medium. Recently Subramanian et al. [15] found experimental evidence that sub-wavelength changes in cellular structure could be detected by analyzing the far-field scattering pattern. It has been concluded that it is possible to differentiate between cells with different fine structure using their far-field scattering patterns.

Because two-photon imaging technique uses a tightly focused beam, the previous plane wave scattering investigations do not offer much direct insight into optical scattering that may occur during two-photon image acquisition. A similar FDTD model with an incident focused beam will aid in development of improved two-photon imaging methods. While such an investigation has not previously been undertaken, the FDTD method has been used in a variety of simulations utilizing a focused beam as an incident source. Such models have been used for investigating the optical trapping force of a focused beam on small particles [16–18] as well as plasmonic generation in nanometallic structures [19, 20].

## 2. Focused Beam FDTD Model

#### 2.1. Scattered-Field Only Formulation

Developed by Yee [9], the FDTD method is a direct solution of Maxwell’s curl equations in the time domain. Maxwell’s curl equations are discretized in space via the central differencing approach and in time by forward differencing. This results in six finite difference equations: one for each component of the electric and magnetic fields. Using these equations in an alternating fashion, the electric and magnetic fields are updated in a leapfrogging manner.

Because the incident field is not included in the FDTD update equations, it must be inserted into the compute space via another method. To minimize memory used to store the field components, we used a scattered-field only formulation [10]. The fields in the compute space are split into incident and scattered field components and conductivity was assumed to be zero, resulting in the following formulation:

The last term in Eq. (1) that depends on the incident field may be redefined as a current source, and the resulting formulation looks identical to the original total-field Yee formulation. At every time step, the incident field is computed at every grid point where *ε _{r}*≠1. Because the source term in Eq. (1) is nonzero only when the grid cell contains a scattering medium, the incident field need not be computed for any grid cells containing the background medium. Because the incident field is computed analytically and does not propagate via the FDTD update equations, it is immune from any phase errors that may occur due to numerical dispersion often associated with FDTD [10].

Because FDTD is by nature a near-field solution, an absorbing boundary condition (ABC) must be applied to the grid edges to reduce reflection from the grid termination. we used the perfectly matched layer (PML) ABC described by Berenger [21], which consists of adding an additional layer of grid cells around the compute space in which each of the field components is split into two other components. In addition, a conductivity term is added to the grid voxels in the PML layer that allows the layer to act as an absorber.

The split field equations allow the conductivity terms to be uniquely defined for each Cartesian direction and the conductivity and PML thickness may be changed to achieve a desired reflectivity from the boundary. This allows a greater flexibility in defining absorbing conditions when compared to other ABCs such as the Mur boundary [22]. These equations can be discretized in the same manner as the original Maxwell’s equations and the resulting 12 update equations can easily be incorporated into the FDTD model. The PML ABC has been extensively used in FDTD simulations involving optical scattering.

#### 2.2. Focused Beam Formulation

The incident beam used in the simulation is a tightly focused Gaussian beam, which is used to closely approximate the incidence conditions that may occur in two photon microscopy. A monochromatic beam propagating in the +*z* direction at a free space wavelength of 800 nm was simulated. The beam waist radius was set to be one wavelength, corresponding to a numerical aperture of 0.5. We note that the focused beam was treated as a continuous wave source. We did not model the temporal properties of the pulsed source, as described in Eq. 4.

While a number of focused beam formulations for use in FDTD models exist [23–25], a formulation developed by Barton and Alexander [25] for a fifth-order corrected Gaussian beam was used in each simulation. This formulation was chosen due to its simplicity, accuracy, and ease of insertion into the model using the scattered-field only method. The equations in Ref. [25] begin by defining a paraxial Gaussian beam that is linearly polarized in the x direction. Because the paraxial beam does not exactly satisfy Maxwell’s equations, additional terms are required in the formulation to reach an acceptably low error level, especially as the beam focal waist is reduced to one wavelength or smaller. These terms add small corrective field components in both the *y* and *z* directions. For a beam waist size of one wavelength (*λ*), the average percent error over all grid locations as calculated in Ref. [25] is approximately 0.01% with a maximum error of less than 1%.

#### 2.3. Cell Configuration and Optical Properties

A three dimensional FDTD code used in a previous study [11] was modified to compute focused beam propagation through regions with complex geometries consisting of cellular components with different refractive indices. An example three cell configuration with each cell containing multiple organelles and a heterogeneous nucleus is shown in Fig. 1. The background is matched to that of cytoplasm (*n*=1.36 [26]) for more accurate modeling of scattering by cells embedded in tissue. The cellular components with higher refractive indices are shown in lighter colors. The dotted outlines show the approximate boundaries of each cuboidal cell. The refractive index values used for the nuclei and mitochondria were *n*=1.40 [26] and *n*=1.38 [27] respectively. The cuboidal cells shown in Figs. 1 and 2 have major diameters of 15 *µ*m and minor diameters of 13 *µ*m. The major and minor diameters of the ellipsoidal nuclei are 6 *µ*m and 5 *µ*m respectively. The organelles are divided evenly by volume into two groups: one ellipsoidal group with major and minor diameters of 1.5 *µ*m and 0.5 *µ*m and another spherical group with diameters of 0.5 *µ*m. The location and orientation of cellular components are chosen via random number generation. The seed for the random number generator can be fixed to allow a particular geometry to be repeated, or it may be changed to allow for original configurations and to test for sensitivity of results to particular cellular configurations.

When included, inhomogeneities within the nucleus (not shown in Fig. 1) consist of randomly placed overlapping spheres of 1 *µ*m diameter. The refractive indices of each sphere are randomly assigned and vary from *n _{n}*-Δ

*n*to

*n*+Δ

_{n}*n*where

*n*=1.4 [26] is the average refractive index of the nucleus and Δ

_{n}*n*is the maximum variation of the refractive index of the inhomogeneities. In all cases Δ

*n*was chosen to be 0.03.

#### 2.4. Simulation Parameters

In order to maintain numerical stability and accuracy using the FDTD method, a grid voxel size of roughly *λ*/20, where *λ* is the shortest wavelength within the scattering medium, is typically used. This means that in order to simulate the optical scattering by a 15 *µ*m biological cell at a wavelength of 588 nm (800 nm light propagating in a *n*=1.36 medium) a problem of one-dimensional size *N*=560 is required. Such a problem would require over 4 GB of memory and roughly 70 processor-hours using the focused beam FDTD code. The three cell simulations require roughly three times the memory and about 500 processor-hours each. Fortunately, FDTD is easily parallelized using the MPI library and can then be run on any number of high performance computing (HPC) systems, such as those at the Texas Advanced Computing Center (TACC) or the Army Research Laboratory Major Shared Resources Center. As an example, a three cell problem can be parallelized over 128 processors and run to completion in about four hours using our program.

Because the source has a turn on time at the first time step, the simulation must be run until the steady state values of the fields are attained. Typically, this requires 3 to 4 times the number of time steps it would take for a signal to propagate from one end of the grid space to another [28]. This results in the computational requirements for FDTD simulations scaling roughly as *N*
^{4}, with *N* as the one dimensional size of the problem grid. Once steady state is reached, the code calculates and records the steady-state optical intensity values along the *x*=0 and *y*=0 planar slices of the 3D space. Additionally, the 3D optical intensity within the focal volume is recorded for accurate calculation of relative two-photon signal strengths.

## 3. Simulations and Results

#### 3.1. Focused Beam Analysis

Planar slices of the steady-state three dimensional optical intensity in the problem space are shown in Fig. 3. The cells in each plot contain a homogeneous nucleus and have 6% organelle density by volume. Fig. 3(a) shows the unscattered incident intensity in the case of a focal depth of approximately 42 *µ*m. Fig. 3(b) shows the total optical intensity for the same incident field. The scattering effect of the cells is evident in the shape of the beam focus and in the increased off-axis signal level in the focal plane. Fig. 3(c) shows the total optical intensity for the case of an incident beam with a focal depth of 22 *µ*m.

To determine the effects of the subcellular components on the PSF, multiple three-cell simulations were run using different cellular configurations. The steady-state optical intensities, examples of which are shown in Fig. 3, were recorded and compared to determine trends in PSF size and two-photon fluorescence excitation. It was also necessary to define metrics for comparing the Gaussian beams for each simulation. The full width half max (FWHM) beamwidths of the squared optical intensity along both the axial and radial directions were compared to those of the incident beam to determine any possible scattering effects on two-photon imaging resolution.

In postprocessing, the axial and lateral profiles of the recorded steady-state Gaussian beam were first fitted to a Gaussian curve of the form seen in Eq. (2). Sample axial and radial beam profiles are shown in Fig. 4(a) and Fig. 4(b), respectively. Note in Fig. 4(a) that the axial profile of the focused beam may vary greatly in half-power beamwidth even when the focal depth is moved by only 10 *µ*m. In Fig. 4(b), it is seen that the radial half-power beamwidth does not not significantly change with focal depth. Increased scattering with depth is seen, however, in the increased optical intensity in the lateral sidelobes of the focused beam. If this trend is extrapolated to greater imaging depths, it is apparent that there will be some maximum imaging depth beyond which this scattered off-axis light will degrade the PSF beyond recovery.

The FWHM was then defined by Eq. (3).

whereσ is taken from Eq. (2). In addition, the recorded steady-state three dimensional optical intensity within the focal volume was used to determine the relative strengths of two-photon fluorescence excitation. The two-photon excitation signal is given by [8]:

We note again that we are only considering the spatial component of the two-photon excitation signal. We are assuming that the temporal dependence of the two-photon excitation signal remains constant with depth. While the two-photon excitation is strictly defined as an integration of the squared optical intensity (*I*
^{2}) over the entire imaged volume, the majority of the two-photon fluorescence will arise from within the beam focal volume. As seen in Eq. (3), the two-photon excitation can therefore be approximated as that excitation that occurs within the focal volume. In this case, the boundaries of the focal volume were defined by the FWHM, as determined in Eq. (3). As seen in the rightmost part of Eq. (4), the intensity within the focal volume was squared, integrated via rectangular quadrature in three dimensions and compared to that of the incident beam to give an indication of the degradation of the induced two-photon fluorescence signal due to scattering.

#### 3.2. Parametric Study

A parametric study was conducted using variations of the multicell geometry seen in Fig. 1. The dotted lines in the figure mark the four focal planes that were used in the study, corresponding to approximate imaging depths of 12, 22, 32 and 42 *µ*m respectively, moving from left to right. For each focal depth, the density of the smaller organelles in each cell was varied between 0% and 10% by volume in 2% increments. The organelle types were distributed evenly by volume between the ellipsoidal scatterers representing mitochondria and the spherical scatterers that represent other organelles such as lysosomes. In addition, the nuclei of the three cells were changed to be either a homogeneous ellipsoid of refractive index of 1.40 or that of the cytoplasm background. This latter refractive index of 1.36 effectively represents the no nucleus case. Placement of the subcellular components in each cell is done via random number generation, but they are constrained to be entirely within each cuboudal cell volume. To investigate the sensitivity of the scattering results to cellular internal configuration, the parametric study was conducted on two different 3-cell geometries, each of which was generated using a different starting random number seed. The total number of simulations initially run was 144. Each three cell simulation required approximately 500 processor-hours.

#### 3.3. Effects of Increasing Focal Depth

It is well-established that two-photon imaging has an in-vivo depth limit of between 300 *µ*m and 1000 *µ*m, depending on a number of factors including tissue type [29]. It is apparent, therefore, that an accurate numerical model of light-tissue interaction should show a loss of two-photon excitation within the focal volume as the depth of the focal plane is increased. In the cases where the cells have no nuclei, this effect is observed as an exponential decrease in two-photon fluorescence excitation over the four increasing focal depths (see Fig. 5(a)).

Fitting the curves in Fig. 5(a) to a simple single scattering attenuation curve of *Ae*
^{-2µsz} yields scattering coefficients of 3.41 *mm*
^{-1}, 6.56 *mm*
^{-1}, 9.99 *mm*
^{-1}, 13.31 *mm*
^{-1} and 16.39 *mm*
^{-1} for organelle volume densities of 2%, 4%, 6%, 8% and 10% respectively. The values were consistent when the placement of the organelles within the cells was varied. These scattering coefficients are also consistent with those reported for dermis or epidermis [30].

#### 3.4. Effects of Volume Organelle Density

In previous numerical studies of light scattering by single inhomogeneous cells it was found that, as expected, small organelles are responsible for a large portion of the large angle far-field scattering [11, 12, 31]. In the focused beam case, this near-isotropic scattering from the smaller organelles was found to cause the predicted two-photon excitation signal to decrease as the volume fraction of organelles in each cell was increased, regardless of focal depth. Intuitively, this decrease in excitation signal with increasing organelle density should become more pronounced as the focal depth is increased, owing to the beam needing to pass through more scattering medium before arriving at the focal plane. This trend is seen in the results seen in Fig. 5(b), which show the predicted two-photon excitation signal as a function of organelle density at each optical depth.

#### 3.5. Effects of Cell Configuration

In the parametric cases involving cells with nuclei, a lensing effect was observed in simulations which placed the focal volume behind or nearly behind the nucleus of a cell. This effect is illustrated in Figs. 3(b) and (c). While this effect was slightly more pronounced in the case of homogeneous nuclei, it remained noticeable even in the case of heterogeneous nuclei as well. The effects of including homogeneous nuclei in the cells are shown in Figs. 6(a) and 6(b). This lensing caused some of the observed trends in two-photon excitation strength to deviate from the results seen in Figs. 5(a) and 5(b). Notably, the nearly linear downward trend in fluorescence excitation within the focal volume seen in Fig. 5(a) is seen at shallow imaging depths in Fig. 6(a), while at greater depths there appears to be an asymptotic trend in the strength of the excitation signal versus focal depth. This suggests that the ellipsoidal nuclei may be acting as a lens and causing a slight recovery of the beam focus at imaging depths greater than 35 *µ*m. This trend may be due to the placement of the nuclei with respect to the axis of the incident beam. When the simulation was re-run using different random number seeds for nuclear placement, it was found that the effect was sensitive to the location of internal cellular structure, particularly as the focal depth was increased. At a focal depth of 22 *µ*m, the fluorescence excitation signal could vary as much as 6% depending on cellular structure. At a depth of 42 *µ*m, the maximum variation increased to 23%. From Fig. 1, it may be inferred that the greater variation at this depth could be due to the nearness of that focal depth to the center of the deepest cell. In certain configurations, this could put the focal plane immediately behind the nucleus of the deepest cell, causing a sharp focusing in these cases and leading to a larger fluctuation in the fluorescence excitation signals from this depth. This lensing effect persisted in the case of cells containing heterogeneous nuclei as well as the homogeneous nuclei case.

In addition, the downward trend in fluorescence signal versus organelle density seen in the no-nucleus case in Fig. 5(b) is drastically reduced by the addition of the homogeneous nuclei. Fig. 6(b) shows that when nuclei are present, the organelle density has relatively little effect on the predicted fluorescence excitation signal when compared their effect when no nuclei are present. This suggests that in this configuration the scattering effects of the nuclei on the two-photon fluorescence signal generated within the focal volume are much greater than those of the organelles.

#### 3.6. Point Spread Function Size

The sizes of the PSF were compared across differing focal depths and varying organelle densities to asses any degradation of spatial resolution. It was found that while there are slight variations of PSF size with increasing depth and organelle density, these variations are typically less than 5% from the incident beam. In some cases the PSF size is actually reduced, improving the theoretical two-photon resolution. Fig. 7(a) shows that the radial resolution of the PSF seems to slightly improve with increasing focal depth, while there is no noticeable trend in PSF size with respect to changes in cell organelle density in Fig. 7(b).

When a homogeneous nucleus is added to each cell, the downward trend in radial PSF width versus focal depth becomes much more pronounced. The drop in PSF size seen in Fig. 8(a) could partially be attributed to the lensing effect of the cell nuclei, particularly at the greatest focal depth. In addition, energy that may just have reached the outer edge of the focal volume in the no-nuclei case could have been scattered further from the focal volume in the homogeneous nuclei case. This increased scattering could then lead to a further reduced PSF size. Fig. 8(b) again shows that PSF size is largely independent of organelle concentration except perhaps for the greatest focal depth case.

#### 3.7. Increasing Number of Cells

Due to computational constraints, the number of cells used in the parametric simulations was kept at three and they were arranged along the z-axis to give the greatest variation in simulated focal depth. However, this configuration seen in Fig. 9(a) does not exactly match what would be encountered by a beam incident on a tissue. Especially in the case of a deeply focused beam, many spatial frequency components of the beam would be affected by cells that are adjacent in the x-y plane to those seen in the simulation. To investigate whether these missing cells would have a significant impact on the results in the study, two additional simulations were run, each with 27 cells arranged in a cube with three cells on each side, as seen in Fig. 9(b). To our knowledge, this is the first such three-dimensional large-scale investigation of optical scattering by multiple heterogeneous biological cells. The one-dimensional size of each 27 cell simulation is *N*=1750, resulting in an approximate memory requirement of 130 GB and a computational requirement of about 5300 processor-hours. The cube of cells is centered on the z-axis and the three cells along that axis are in the exact configuration as those using one of the random seeds in the parametric study. This allows a more direct comparison to be made between the three cell and 27 cell cases. A cell organelle density of 6% by volume was chosen and the beam focus was placed at depths of both 12 *µ*m and 42 *µ*m, corresponding to the shallowest and deepest of the focal planes shown in Fig. 1.

Two-dimensional slices of the steady-state optical intensity data for both the 3 cell and 27 cell cases are shown in Fig. 10. The black dotted box in Fig. 10(b) shows the location of the 3 cell configuration within the 27 cell cube. Figs. 11(a) and 11(b) show those same planar intensity slices in the context of the three-dimensional scattering geometry. A similar analysis of this intensity data was used to determine axial and lateral beam profiles for each case. As might be expected, Figs. 12(a) and 12(b) show that at a focal depth of 12 *µ*m the axial and radial beam profiles are nearly identical, suggesting that off-axis scatterers do not play a significant role in determining focal characteristics for shallow imaging depths. As may also be expected, at deeper focal depths the off-axis cells have a larger role in scattering the incident light. The three cell case resulted in radial PSF beamwidths of 107.1% and 72.6% of the incident width, respectively, for the 12 *µ*m and 42 *µ*m focal depths. The 27 cell case resulted in radial beamwidths of 107.1% and 66.0% of the incident width, showing good agreement with the three cell case and suggesting that additional scattering by off-axis cells does not cause an appreciable degradation in the size of the two-photon imaging PSF as seen in Fig. 12(b). This is consistant with reported measurements of two photon PSF trends with depth which state that radial PSF width is not noticably increased by scattering within the imaging medium within the first 100 *µ*m of imaging depth [4, 5].

In the three cell case, the predicted two-photon excitation signals for focal depths of 12 *µ*m and 42 *µ*m were 71.6% and 18.0% of the incident signal, respectively. Those same signals in the 27 cell case were 71.6% and 9.8%., suggesting that off-axis scatterers contribute to the reduction of two-photon excitation signal in the focal volume. This is again consistant with reported measurements of two-photon fluorescence signal versus depth [5].

These results are also consistent with the probabilistic model of photon scattering, since it is assumed that only ballistic photons will arrive in the focal volume while scattered photons will typically be deflected at an angle large enough to miss the focal volume entirely [5]. This will result in the two-photon PSF remaining relatively unchanged in size, while the fluorescence signal generated within the focal volume will decrease with depth or due to any increase in scattering.

## 4. Conclusion

The FDTD method has been applied to focused Gaussian optical scattering from multiple heterogeneous cells. To our knowledge this is the first three dimensional computational investigation of a focused beam interacting with multiple biological cells using FDTD. A parametric study was performed whereby the simulated cells were varied by organelle density, nuclear type and arrangement of internal cellular structure and the beam focus depth was varied within the cluster of cells. Results showed that scattering is responsible for a reduction of two-photon excitation signal due to reduced optical intensity within the focal volume. Of the organelle types investigated, it appears that the cell nuclei are responsible for the greatest scattering of the focused beam in the configurations studied. No significant degradation of two-photon imaging resolution was observed to occur within the first 40 *µ*m of imaging depth. The results also demonstrate the need for an accurate method to measure the refractive indices of subcellular components to facilitate more accurate modeling results.

## Acknowledgments

The authors would like to acknowledge funding from the National Science Foundation under Grant CBET/0737731, the American Heart Association under Grant 0735136N, the Dana Foundation and the Consortium Research Fellows Program with the U.S. Air Force Research Laboratory/Human Effectiveness Directorate. The authors also acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing high performance computing and visualization resources that have contributed to the research results reported within this paper. URL: http://www.tacc.utexas.edu. Additional computing resources were provided by collaborators with the Directed Energy Bioeffects Division at the U.S. Air Force Research Laboratory at Brooks City-Base in San Antonio, Texas and the Army Research Laboratory Major Shared Resources Center.

## References and links

**1. **W. Denk, J. Strickler, and W. Webb, “2-Photon Laser Scanning Fluorescence Microscopy,” Science **248(4951)**, 73–76 (1990). [CrossRef]

**2. **W. Denk and K. Svoboda, “Photon upmanship: Why multiphoton imaging is more than a gimmick,” Neuron **18(3)**, 351–357 (1997). [CrossRef]

**3. **X. Deng and M. Gu, “Penetration depth of single-, two-, and three-photon fluorescence microscopic imaging through human cortex structures: Monte Carlo simulation,” Appl. Opt. **42(16)**, 3321–3329 (2003). [CrossRef]

**4. **C.-Y. Dong, K. Koenig, and P. So, “Characterizing point spread functions of two-photon fluorescence microscopy in turbid medium,” J. Biomed. Opt. **8(3)**, 450–459 (2003). [CrossRef]

**5. **A. K. Dunn, V. P. Wallace, M. Coleno, M. W. Berns, and B. J. Tromberg, “Influence of optical properties on two-photon fluorescence imaging in turbid samples,” Appl. Opt. **39(7)**, 1194–1201 (2000). [CrossRef]

**6. **X. Gan and M. Gu, “Effective point-spread function for fast image modeling and processing in microscopic imaging through turbid media,” Opt. Lett. **24(11)**, 741–743 (1999). [CrossRef]

**7. **M. Oheim, E. Beaurepaire, E. Chaigneau, J. Mertz, and S. Charpak, “Two-photon microscopy in brain tissue: parameters influencing the imaging depth,” J. Neurosci. Methods **112(2)**, 205–205 (2001). [CrossRef]

**8. **P. Theer and W. Denk, “On the fundamental imaging-depth limit in two-photon microscopy,” J. Opt. Soc. Am. A **23(12)**, 3139–3149 (2006). [CrossRef]

**9. **K. Yee, “Numerical solution of inital boundary value problems involving maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. **14(3)**, 302–307 (1966).

**10. **A. Taflove, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, 3rd ed. (Artech House Inc., Norwood MA, 2005).

**11. **A. Dunn and R. Richards-Kortum, “Three-dimensional computation of light scattering from cells,” IEEE J. Sel. Top. Quantum Electron. **2(4)**, 898–905 (1996).

**12. **C. Liu and C. E. Capjack, “Effects of cellular fine structure on scattered light pattern,” IEEE Trans. Nanobiosci. **5(2)**, 76–82 (2006). [CrossRef]

**13. **J. Q. Lu, P. Yang, and X. H. Hu, “Simulations of light scattering from a biconcave red blood cell using the finite-difference time-domain method,” J. Biomed. Opt. **10**, 024,022 (2005). [CrossRef]

**14. **X. Li, A. Taflove, and V. Backman, “Recent progress in exact and reduced-order modeling of light-scattering properties of complex structures,” IEEE J. Sel. Top. Quantum Electron. **11(4)**, 759–765 (2005).

**15. **H. Subramanian, P. Pradhan, Y. Liu, I. R. Capoglu, X. Li, J. D. Rogers, A. Heifetz, D. Kunte, H. K. Roy, A. Taflove, and V. Backman, “Optical methodology for detecting histologically unapparent nanoscale conse-quences of genetic alterations in biological cells,” Proc. National Acad. Sci. **105(51)**, 20,118–20,123 (2008).

**16. **R. Gauthier, “Computation of the optical trapping force using an FDTD based technique,” Opt. Express **13(10)**, 3707–3718 (2005). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-10-3707. [CrossRef]

**17. **T. A. Nieminen, H. Rubinsztein-Dunlop, N. R. Heckenberg, and A. I. Bishop, “Numerical Modelling of Optical Trapping,” Comp. Phys. Comm. **142**, 468–471 (2001). [CrossRef]

**18. **W. Sun, S. Pan, and Y. Jiang, “Computation of the optical trapping force on small particles illuminated with a focused light beam using a FDTD method,” J. Mod. Opt. **53(18)**, 2691–2700 (2006). [CrossRef]

**19. **K. Choi, H. Kim, Y. Lim, S. Kim, and B. Lee, “Analytic design and visualization of multiple surface plasmon resonance excitation using angular spectrum decomposition for a Gaussian input beam,” Opt. Express **13(22)**, 8866–8874 (2005). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-22-8866. [CrossRef]

**20. **R.W. Ziolkowski, “FDTD modeling of Gaussian beam interactions with metallic and dielectric nano-structures,” in *Proc. 2004 URSI International Symposium on Electromagnetic Theory*, pp. 27–29 (2004).

**21. **J.-P. Berenger, “Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves,” J. Comput. Phys. **127(2)**, 363–379 (1996). [CrossRef]

**22. **G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic field equations,” IEEE Trans. Electromagn. Compat. **23(4)**, 377–382 (1981). [CrossRef]

**23. **I. R. Capoglu, A. Taflove, and V. Backman, “Generation of an incident focused light pulse in FDTD,” Opt. Express **16(23)**, 19,208–19,220 (2008). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-23-19208.

**24. **W. Challener, I. Sendur, and C. Peng, “Scattered field formulation of finite difference time domain for a focused light beam in dense media with lossy materials,” Opt. Express **11(23)**, 3160–3170 (2003). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-11-23-3160 [CrossRef]

**25. **J. P. Barton and D. R. Alexander, “Fifth-order corrected electromagnetic field components for a fundamental Gaussian beam,” J. Appl. Phys. **66**, 2800 (1989). [CrossRef]

**26. **A. Brunsting and P. F. Mullaney, “Differential Light Scattering from Spherical Mammalian Cells,” Biophys. J. **14(6)**, 439 (1974). [CrossRef]

**27. **H. Liu, B. Beauvoit, M. Kimura, and B. Chance, “Dependence of tissue optical properties on solute-induced changes in refractive index and osmolarity,” J. Biomed. Opt. **1(2)**, 200–211 (1996). [CrossRef]

**28. **A. Taflove and M. Brodwin, “Numerical Solution of Steady-State Electromagnetic Scattering Problems Using the Time-Dependent Maxwell’s Equations,” IEEE Trans. Microwave Theory Tech. **23(8)**, 623–630 (1975). [CrossRef]

**29. **M. Oheim, D. J. Michael, M. Geisbauer, D. Madsen, and R. H. Chow, “Principles of two-photon excitation fluorescence microscopy and other nonlinear imaging approaches,” Adv. Drug Del. Rev. **58(7)**, 788–808 (2006). [CrossRef]

**30. **A. J. Welch and M. J. C. van Gemert, eds., *Optical-Thermal Response of Laser-Irradiated Tissue*, 1st ed. (Plenum Press, 1995).

**31. **R. Drezek, A. Dunn, and R. Richards-Kortum, “Light scattering from cells: finite-difference time-domain simulations and goniometric measurements,” Appl. Opt. **38(16)**, 3651–3661 (1998).