A Pulsed Finite-difference Time-domain (fdtd) Method for Calculating Light Scattering from Biological Cells over Broad Wavelength Ranges

We combine the finite-difference time-domain method with pulse response techniques in order to calculate the light scattering properties of biological cells over a range of wavelengths simultaneously. The method we describe can be used to compute the scattering patterns of cells containing multiple heterogeneous organelles, providing greater geometric flexibility than Mie theory solutions. Using a desktop computer, we calculate the scattering patterns for common homogeneous models of biological cells and also for more complex representations of cellular morphology. We find that the geometry chosen significantly impacts scattering properties, emphasizing the need for careful consideration of appropriate theoretical models of cellular scattering and for accurate microscopic determination of optical properties. In vivo spectrophotometric evaluation of neoplastic and non-neoplastic skin pigmented lesions. Spectroscopic diagnosis of bladder cancer with elastic light scattering, " Lasers Surg. Elastic scattering spectroscopy as a diagnostic for differentiating pathologies in the gastrointestinal tract: preliminary testing, " J. Identification of colonic dysplasia and neoplasia by diffuse reflectance spectroscopy and pattern recognition techniques, " Appl. reflectance spectroscopy of human adenomatous colon polyps in vivo, " Appl. Predictions and measurements of scattering and absorption over broad wavelength ranges in tissue phantoms, " Appl. Observation of periodic fine structure in reflectance from biological tissue: a new technique for measuring nuclear size distribution, " Opt. Lett. 80, 627-630. Polarized light scattering spectroscopy for quantitative measurement of epithelial structures in situ, " IEEE J.Reflectance spectroscopy with polarized light: is it sensitive to cellular and nuclear morphology, " Opt. Mechanisms of light scattering from biological cells relevant to noninvasive optical-tissue diagnostics, " Appl. Light scattering from coated spheres: model for biological cells, " Appl. A transmitting boundary for transient wave analysis, " Sci. Sin. Light scattering from cells: finite-difference time-domain simulations and goniometric measurements, " Appl. Improvements to the finite-difference time-domain method for calculating the radar cross section of a perfectly conducting target, " IEEE Trans. Microwave Theory Tech. 38, 919-927 (1990). 21. C. Britt, " Solution of electromagnetic scattering problems using time domain techniques, " IEEE Trans.Light scattering from biological cells: dependence of backscatter radiation on membrane thickness and refractive index," Appl. Fractal texture features based on optical density surface area: use in image analysis of cervical cells, " Analyt.


Introduction
Recent work has suggested that elastic light scattering spectroscopy can provide a valuable, non-invasive means to quantitatively probe tissue morphology.A number of studies have demonstrated that elastic scattering can be used to differentiate normal and diseased tissues in several organ sites including the skin, bladder, and colon [1][2][3][4][5].More recently, several investigators have reported methods to extract quantitative features related to tissue morphology from reflectance spectra.Mourant et al. reported that the size of scatterers could be determined from elastic scattering spectra using Mie theory [6].This work involved tissue phantoms comprised of polystyrene spheres and Intralipid.Perelman et al. described a technique for obtaining nuclear size distributions by Fourier transformation of the scattering spectra [7].Zonios et al. reported a method to extract effective scatterer sizes from colon tissue spectra based on a Mie theory model [5].In order to isolate the scattering signal due to cell nuclei from diffuse background scattering, which is typically orders of magnitude larger, Backman et al. [8] and Sokolov et al. [9] proposed techniques based on polarized light scattering spectroscopy.These polarization-based techniques assessed size and refractive indices of cell nuclei using Mie theory models.The use of polarization provided a means to distinguish between single and multiply scattered light and reduced the effects of hemoglobin absorption.All of the studies noted above developed models based on Mie theory which could adequately describe experimental data, offering a potential approach to extracting quantitative information from the scattering spectra.
Although a number of studies have shown promising results using Mie theory models, the ability to match a model with experimental data does not necessarily indicate the model is the best model or even an appropriate one.There is not convincing experimental evidence that justifies the use of Mie theory as an accurate model of cellular and/or nuclear scattering.Mie theory is commonly applied to model scattering from cells largely because it is the only simple analytical model available.In fact, experimental evidence at times suggests a Mie theory model of a biological cell may not be appropriate.For instance, work by Mourant et al. found that the scatterer sizes in biological cells ranged from 0.4 to 2.0 µm, a size consistent with organelles smaller than a cell's nucleus [10].Moreover, the same study showed that the nucleus contributed mostly to low angle scattering while smaller organelles contributed to high angle scattering, which may be particularly pertinent to reflectance measurements acquired using fiber-based probes with small source-detector separations.Other studies have also produced results which indicate that Mie theory may not always be a sufficient model.For instance, McGann et al. clearly demonstrated that forward scattering of lymphocytes varied inversely with cell volume.This would not be expected from a Mie theory model of scattering and emphasizes that cellular light scattering cannot always be adequately described using the simplest conceivable geometric model [11].Additional research has described significant changes in backscattering obtained when a coated sphere (cell with membrane) [12] or concentric sphere (nucleated cell) [13] model was used rather than a single sphere model, further emphasizing the importance of a prudent approach to choosing a cell model since the scattering predicted is a strong function of the geometry chosen.
Based on the available experimental evidence, it seems likely that there will be cells which cannot be adequately modeled as a homogeneous sphere.For such cases, a more flexible alternative than Mie theory is required.We provide one possible solution in this paper.The method we describe can be used to calculate light scattering from cells of arbitrarily complicated shape and dielectric structure.Our approach is based on a standard finite-difference time-domain (FDTD) technique [14].In the past, the primary limitation of the FDTD technique has been the need for an individual FDTD run for each frequency of interest, making the technique impractical for calculations at more than a few frequencies.In this paper, we modify the standard approach incorporating pulse response techniques which provide a means to calculate light scattering properties over a broad frequency range using only a single FDTD run.
Because the computational time requirements of the proposed method are extremely modest, requiring little more time than equivalent Mie theory calculations, the method can be efficiently used to compute light scattering for any desired model of a cell, from models as simple as a single sphere or coated sphere to cell models containing numerous organelles and intricate dielectric structure.The purpose of this paper is two-fold: first, to provide a detailed description of the practical implementation of the method, and second, to use the method to compute the scattering properties of cells modeled using progressively more complex geometries and to demonstrate that the cell model used significantly impacts the scattering properties.

Methods
The FDTD algorithm provides a means to numerically solve Maxwell's equations in the time domain.The method is briefly reviewed here.A more detailed explanation is found in [14].The algorithm begins by discretizing Maxwell's curl equations in space and time, resulting in a set of explicit finite-difference equations.The finite-difference equations are stepped in time, and the electric and magnetic field components at each grid point are alternately updated.To prevent artificial reflections along the edges of the grid, an appropriate boundary condition must be employed.In this paper, the Liao boundary condition was used [15].For applications requiring a larger dynamic range, the PML boundary condition may be applied to provide at least 80 dB additional suppression [16].
To yield accurate results for a given wavelength, the grid spacing used must be less than wavelength, typically λ/10 or smaller.At each grid point, the permittivity and conductivity of the medium is specified.The cell is constructed by assigning permittivity values to each cell component.A range of values may be assigned to a particular component if that component is inhomogeneous.Details regarding the application of the FDTD method specifically to biological cells, including refractive index values for specific cellular constituents, may be found in [17][18][19].
In a standard FDTD model, a plane wave is propagated throughout the grid until a steady state solution is obtained.At this point, the electric and magnetic field values are known on the entire grid which lies in the near field.To compute the scattering pattern in the far field, a near field to far field transformation is required.For two-dimensional FDTD, this is accomplished by weighting the near field data with the free space Green's function, written in terms of the zero-order Hankel function, H 0 (1) (r), and integrating over a surface completely encompassing the cell [14].Using this conventional approach, it is necessary to complete a separate FDTD run for each frequency of interest.
A number of researchers have proposed more efficient approaches to obtaining a frequency response using the FDTD method [20,21].These approaches have been applied to calculate the frequency response of such objects as perfectly conducting cylinders, cones, spheres, and boxes [20,21].In this paper, we apply similar methods to calculate the frequency response over a range of interest for biological cells.The fundamental idea behind frequency response FDTD is to provide excitation in the form of a time-limited pulse and to obtain a frequency response by Fourier transformation of the field component waveforms as described below [20].
The pulsed FDTD technique provides a means to obtain a frequency response in only one FDTD run with approximately the same computational requirements of a single frequency run.The grid dimensions and incident pulse are determined based on the lowest wavelength of interest.Either a Gaussian or raised cosine pulse can be used as the excitation source for pulsed FDTD; in this work, a Gaussian pulse was employed.To obtain the frequency response, the magnitude and phase of the time domain waveforms are determined using a discrete Fourier transform (DFT).The discrete Fourier transform is expressed as shown in Eq. ( 1): where g(n∆t) is the field value, n is the time step, N is the length of the DFT, ∆f is the frequency resolution, k is the frequency index, and NF is the number of frequencies [20].The summation in Eq. 1 is updated at each time step for each field component.At the end of the simulation, the frequency values are normalized by the DFT of the incident pulse.The problem converges when the applied pulse has decayed to zero for all grid positions.
It should be pointed out that for the two-dimensional case, the memory required for continuous wave (single frequency) FDTD increases linearly with the number of grid elements.Given a grid containing 250x250 elements, the memory required for a single frequency calculation is approximately 500 Kb to store the necessary field values and 0.18 Kb to store far field intensity information sampled every 2 degrees.To perform the same simulation using 100 frequencies, 500 Kb is still required to store field values, but now 100 times more memory, or 18 Kb, is required to store the far field intensities.Relative to the overall problem size, the additional memory burden imposed to obtain spectral data rather than single frequency data is very minor.All simulations described in this paper were performed using a two-dimensional FDTD code.The 2D code is advantageous in that memory requirements are modest (<1 Mb of memory) and calculations are rapid, requiring approximately two minutes on a 500 MHz desktop computer.The techniques described are equally applicable to three-dimensional simulations.Since only sparse references to pulsed FDTD are found in the electromagnetic literature, a series of simulations were run to verify that the method produced results in agreement with Mie theory for a series of non-conducting objects.The first test geometry consisted of a 12 µm grid containing a 4 µm circular object (m = 1.02).The scattering diagrams for wavelengths ranging from 1 µm to 2 µm with a 5 nm increment were computed.The use of the terms scattering diagram, a, n, x, and m follow the definitions of van de Hulst [22].The term "scattering diagram" refers to a plot of the scattered intensity as a function of angle.

Code Verification
When the scattering diagram is normalized so that a value of one is obtained when integrated over 4π steradians, the curve is referred to as a phase function.Throughout this paper, the quantity m refers to relative refractive index, a refers to the scatterer radius, λ = λ 0 /n where n is the index of refraction of the medium surrounding the object and λ 0 is the wavelength in vacuum, and x = 2πa/λ.FDTD data were compared to an analytical solution computed using Mie theory.Results are shown in Figure 1.The FDTD curves agree closely with theoretical predictions.The discrepancy between Mie theory and FDTD data at the highest angles is due to imperfect boundary conditions.The two curves shown in Figure 1 demonstrate that in some cases there is a precise match with theoretical predictions at all angles while in other instances imperfections in the boundary conditions can influence the high angle data (>165°).The top curve in Figure 1 shows worst case data; the bottom curve shows best case data.It should be pointed out that high angle discrepancies between FDTD data and theory are most prominent when the dynamic range of the problem is large.Curves are normalized to scattered intensity at 0°.
To provide a more extensive test of the code against Mie theory, the FDTD code was used to replicate curves shown in van de Hulst [22].The simulation used a homogeneous cylinder geometry (m=1.5) with x ranging from 1.2 to 2.4.Tabulated values for particular points along the set of curves {x=1.2, x=1.6, x=2.0, x=2.4} were calculated by Rayleigh and are tabulated in [22].FDTD results are shown in Figure 2. The left half of the figure shows four individual scattering diagrams.The right half of the figure shows the complete set of results.The vertical axis of each scattering diagrams is the scattered intensity, or using the notation of van de Hulst, the quantity (2/π)|T(θ)| 2 [22].There were not high angle discrepancies in this simulation.

Comparison of Cell Geometries
To demonstrate some of the potential applications of pulsed FDTD to modeling scattering from biological cells, a number of simulations were performed, each over a broad wavelength range.The first set of simulations was designed to demonstrate how the different organelles which constitute a cell contribute to scattering.Four cases were considered.First, the crosssection of the nucleus of a cell was modeled as a homogeneous circle.Second, the cytoplasm of the cell was modeled as a homogenous circle.These homogeneous cases are the geometries most often used to generate data concerning cellular scattering since the scattering intensities are obtainable via Mie theory.After considering the first two cases, the cell was next modeled using a concentric circle geometry consisting of a homogeneous circle representing the cytoplasm and a second interior circle representing the nucleus.Several investigators have employed a concentric sphere model of cellular scattering since analytical solutions for this geometry are available in the literature [12,23].Finally, organelles of various sizes and indices of refraction were added into the concentric circle model (cytoplasm plus nucleus) to demonstrate the effect of heterogeneities in refractive index structure on scattering properties.This set of simulations, in which cellular geometries were defined in a progressively complex manner, is useful in examining how various portions of the cell contribute to scattering.
For each of the four geometries described, the light scattering intensity functions (scattering diagrams) were calculated over a range of wavelengths.To facilitate comparisons between the four cases, the cell cytoplasm when present had a diameter of 8 µm, and the Wavelengths spanned from 600 nm to 1000 nm with a 5 nm increment.Results are shown in Figure 3.All of the homogeneous cell models produced scattering diagrams with clear interference peaks.Differences in the spacing of the peaks based on the size of the scatterers are clearly seen.Spreading of the peaks as the scatterer size relative to the wavelength increases is evident.Although at quick glance the four data sets may appear similar, the details of each case are quite different.It is particularly interesting to note how the introduction of heterogeneities in the form of small organelles impacts scattering.The addition of cytoplasmic organelles begins to obscure the interference peaks clearly visible in the simulations using homogeneous geometries.The effects of the heterogeneities are most noticeable at angles over ninety degrees, partially because the scattered intensity values in this region are five to six orders of magnitude smaller than the scattered intensity values at low angles.The trends in high angle scattered intensity as a function of wavelength using the different cellular models are also significant, demonstrating that the wavelengths and scatterer sizes relevant to optical diagnostics often lie in a scattering regime in which high angle scatter is particularly sensitive to scatterer size.These types of observations underscore the necessity of careful consideration of which portions of a cell are responsible for observed scattering before choosing a simple geometric model and then using the model to create theoretical scattering curves to be fit to experimental data.

Heterogeneous Cell Geometries: Normal versus Pre-cancerous Cells
To consider how more complicated descriptions of cellular morphology influence light scattering properties, two cells containing multiple sizes and shapes of organelles and heterogeneous nuclei were considered.In the first cell, the morphology was defined using histological features of normal cervical cells; in the second cell, the morphology was defined based on the features of cervical cells staged as high grade dysplasia [24][25][26].In order to emphasize differences due to the internal contents of the two cells rather than the overall cell size and shape, both cells were the same size (9 µm diameter) and circular in shape.The most significant differences between the dysplastic cell relative to the normal cell included increased nuclear size and nuclear-to-cytoplasmic ratio (normal 0.2, dysplastic 0.67), asymmetrical nuclear shape, increased DNA content, and hyperchromatic nucleus with areas of coarse chromatin clumping and clearing.In the normal cell, nuclear index variations were uniformly distributed between ∆n=±0.02 about the mean nuclear index, n=1.40, at spatial frequencies ranging from 10-30 µm -1 providing a fine, heterogeneous chromatin structure.In the dysplastic cell, nuclear index variations were distributed between ∆n=±0.04 about the mean nuclear index, n=1.42, at spatial frequencies ranging from 3-30 µm -1 providing a coarser, more heterogeneous chromatin structure.Both normal and dysplastic cells contained several hundred organelles (radii ranging from 50 nm to 0.5 µm; n=1.38 to 1.40) randomly distributed throughout the cytoplasm.Wavelengths spanned from 600 nm to 1000 nm with a 5 nm increment.
Results are shown in Figure 4. Elevated scattering is evident in the dysplastic cell.The increased scattering at small angles is due to the larger nucleus.The increased scattering at higher angles is due to alterations in chromatin structure, resulting in increased heterogeneity of refractive index variations.Since the dysplastic cell contains a large heterogeneous nucleus containing an assortment of scatterer sizes and refractive indexes, distinct interference peaks are not present.In the normal cell, although heterogeneities are present in the cellular structure, they are not significant enough to disrupt the peaks resulting from the cytoplasm and nuclear boundaries.
The bottom half of Figure 4 displays the integrated scattered intensity as a function of wavelength for three angular regions: 0-20°, 80-100°, and 160-180°.The curves show that the integrated intensity is a function of both angle and cellular structure.In this simulation, changes in the wavelength dependence of the scattering between the normal and dysplastic cell are particularly evident at high angles.To develop optimized optical techniques which can discriminate between normal and dysplastic tissue based on differences in the wavelength dependence of cellular scattering, it is important to investigate which angular regions offer the highest potential for differential diagnosis.Then optical probes with delivery and collection geometries designed to preferentially sample particular angular ranges can be developed.

Discussion and Conclusions
The scattering properties of cells modeled using homogeneous and heterogeneous descriptions of morphology are clearly, and markedly, different.However, the goal of this work is not to suggest that either a homogeneous or heterogeneous model of a biological cell is more advantageous than the other, but rather to offer a novel approach that can be applied in an equally straightforward manner regardless of the proposed geometry.Using such an approach, it is possible to consider those representations of cells which experimental evidence suggests are most sensible, rather than imposing artificial limits on potential geometries based on the availability of simple analytical solutions.Although the FDTD method offers tremendous flexibility and is widely employed for electromagnetic modeling applications in areas such as antenna design, microelectronics, and the defense industry [14], it has not found extensive use in the field of biomedical optics.Perhaps this is partly due to the computational requirements of the method which until the past several years limited simulations of biological cells at optical frequencies to large supercomputers.Today the same problems can be solved in both two and three dimensions on desktop computers.The FORTRAN program used for the simulations described in this paper requires approximately two minutes for a two-dimensional solution and ninety minutes for a three-dimensional solution using a 500 MHz workstation.Using double precision variables, memory requirements for a 10 µm cell are approximately 1 Mb in two dimensions and 60 Mb in three dimensions.
Until the last few years, an additional limitation of FDTD calculations was the lack of sufficient boundary conditions to adequately suppress reflections of outgoing waves at grid boundaries.Imperfect boundary conditions, in certain cases, can limit the accuracy of results  [16], have almost totally eliminated the occurrence of incompletely terminated outgoing waves.The new boundary conditions are based on surrounding the FDTD problem space by a "perfectly matched" absorbing material which perfectly transmits waves of all polarizations, all frequencies, and all angles of incidence [27].The combination of significant improvements in absorbing boundary conditions coupled with increased available computational speed and memory makes the FDTD technique more practical than ever before.
In previous work, we demonstrated that the choice of cell geometry has a profound impact on scattering at a particular frequency [19].In particular, when the cell is treated as a structure containing a continuum of refractive index fluctuations rather than a structure composed of several homogeneous elements, the obtained scattering data are strikingly altered.Although the potential to use FDTD as an alternative to Mie theory has existed for several years, without a method to efficiently perform FDTD computations at a large number of frequencies, it was not realistic to extend the use of FDTD to the problem of extracting morphologic parameters from reflectance data.By incorporating pulse response techniques into the FDTD algorithm, the code can be used to perform the types of calculations required in the process of inverting reflectance data as described in [5,[7][8][9].
Because it is possible to obtain such varied pictures of cellular scattering based on the assumptions regarding cell geometry, it is crucial to pursue the fundamental experimental work necessary to elucidate the sub-cellular sources of scattering in tissue.Can epithelial tissue really be treated as a collection of spheres in which all scatters except for the nucleus are neglected, or do the cell membranes and other organelles in between significantly matter?Typically, when observing images of tissue, one views histological micrographs in which the visibility of the nuclei is artificially enhanced by exposing the nuclei to absorbing stains.Observing the same tissue using a reflectance confocal microscope or other scattering-based imaging method, the picture obtained can be much more complicated.Depending on the type of the cell, the nucleus may or may not be a significant scatterer.Even if the nucleus is the dominant scatterer, is it the nucleus taken as a homogeneous entity or is it the size and refractive index of the chromatin particles within it that are really dictating the observed scattering properties?Although computational studies can provide results given a particular set of simulations parameters, only experimental work can determine what the truly significant scatterers are for a particular tissue and how the scatterers change during the course of disease progression.As further experimental work is conducted to provide a clearer understanding of cellular scattering, it will become increasingly feasible to use methods such as pulsed FDTD to perform the calculations necessary to accurately extract morphologic data from tissue spectra.
It is important to emphasize that the results presented in this paper were obtained using a two-dimensional FDTD model.Implementing the pulse response techniques in a twodimensional model is a much simpler problem than the three-dimensional case.Given the high level of interest throughout the past year in developing methods to extract morphological parameters from reflectance data, we believe it important to point out the potential benefits afforded by pulsed FDTD approach in a timely fashion.In future work, the technique will be extended to the three-dimensional case.Two-dimensional FDTD modeling of cells has been explored previously and work comparing 2D and 3D solutions of the problem indicated qualitatively similar trends [28].The influence of parameters such as increasing nuclear to cytoplasmic ratio and adding different sizes and shapes of organelles produced similar results in the two and three dimensional cases.To incorporate absolute scattered intensities into a model which would ultimately be applied to the analysis of measured reflectance data, a threedimensional pulsed FDTD approach should be used.A three-dimensional implementation of pulsed FDTD offers a new approach to modeling cellular scattering for applications which attempt to assess morphology via reflectance measurements.Pulsed FDTD has other potential applications as well such as generating improved phase functions for Monte Carlo simulations or developing novel analysis methods for flow cytometry.

Fig 2 .
Fig 2. Validation of the pulse response FDTD code.Calculated scattering diagrams for a series of infinite cylinders (m=1.50, x=1.2 to 2.4).The left graph shows curves obtained by plotting the intensity data corresponding to particular size parameters.The right image displays the calculated scattering over a range of size parameters.The color scale corresponds to the log of the scattered intensity.Each curve in the graph on the left corresponds to a horizontal line through the image on the right.

4 (
C) 2000 OSA 27 March 2000 / Vol. 6, No. 7 / OPTICS EXPRESS 152nucleus had a diameter of 4 µm.Refractive index values for the cytoplasm and the nucleus were 1.37 and 1.40, respectively.Organelle refractive indices ranged from 1.38 to 1.42, and organelle sizes ranged from 100 nm to 1 µm.Approximately 25% of the available space within the cell (space not already occupied by the nucleus) was filled with organelles.

Fig 3 .
Fig 3. Four models of cellular scattering: (1) nucleus only, (2) cytoplasm only, (3) nucleus and cytoplasm, and (4) nucleus and cytoplasm containing organelles.The color scale corresponds to the log of the scattered intensity.

Fig 4 .
Fig 4. Top: Scattering from normal (left) and dysplastic (right) cervical cell.Note elevated scattering in dysplastic cell.The increased scattering at small angles is due to a larger nucleus to cytoplasm ratio.The increased scattering at high angles is due to alterations in chromatin structure, resulting in increased heterogeneity in nuclear refractive index.The color scale corresponds to the log of the scattered intensity.Bottom: Integrated scattered intensities over three angular ranges (0-20°, 80-100°, and 160-180°) for normal (left) and dysplastic (right) cervical cells (C) 2000 OSA 27 March 2000 / Vol. 6, No. 7 / OPTICS EXPRESS 155 at scattering angles approaching 180 degrees.However, boundary conditions have improved dramatically over the past five years.The newest boundary conditions, all variants of the perfectly matched layer (PML) solution developed by Berenger