1.

Introduction

Understanding the nature of light interaction with tissue is fundamental in aiding the development of optical diagnostic techniques capable of fast, noninvasive evaluation of tissue pathology. A number of these modalities, such as diffuse optical spectroscopy¹ and elastic scattering spectroscopy,² depend on correlating the scattering properties of tissue to its physiological state. Scattering spectroscopy techniques have been developed that can determine the degree of tissue dysplasia in an effort to diagnose various forms of cancer.^{3, 4, 5} Many of these techniques would benefit greatly from a large set of training data that spans a large parameter space of cell variation. While such a data set may be difficult to generate experimentally, a computational method would enable the user to independently control all physical parameters of the tissue, such as cell size, nuclear size, or organelle density.

A direct analytical solution to optical scattering problems involving cells with any realistic level of complexity is virtually impossible. Analytical solutions to scattering problems such as Mie theory⁶ and variations of it involving coated spheres⁷ are available, but they are often very limited in the situations they may be appropriately applied to. While spectroscopic techniques based on Mie theory have demonstrated some ability to determine scatterer sizes from elastic scattering data,^{4, 8} all such approximations require a significant simplification of the cell geometry and refractive indices of cell components.

With the ever-increasing capabilities of high-performance computing resources, it is possible to solve complex scattering problems by treating Maxwell’s equations numerically using computational techniques such as the finite-difference time-domain (FDTD) method. The FDTD method has been used extensively to investigate optical scattering from single heterogeneous cells^{9, 10, 11, 12, 13, 14, 15, 16, 17, 18} as well as scattering from multiple homogeneous cells.^{19, 20, 21} It has recently been demonstrated that optical scattering problems involving up to 27 complex cells are feasible by running a parallelized FDTD code on high-performance computing clusters.²²

Such a large problem requires access to significant computational hardware and a high number of CPU hours. With the increasing number of open scientific computing resources, such as the National Science Foundation’s Teragrid,²³ access to such computing power is available, but resources are still limited. Therefore, it is of interest to reduce the number of clock cycles required to run a simulation without significantly degrading the simulation results.

A linear superposition method that is a modified version of the near-to-far-field transform is presented in this paper. It is evaluated for its effectiveness in estimating a large multicell far-field scattering pattern from the superimposed far-field patterns of single cells or small numbers of complex cells. While it has previously been demonstrated as a method to estimate scattering of small numbers of red blood cells from that of a single cell,¹⁹ it is unclear whether such a technique will be as effective on a complex group of scatterers such as heterogeneous cells. The objective of this paper is to clarify the requirements and limitations of the superposition method in estimating the simulated farfield scattering from multiple biological cells from simulations of a reduced number of cells.

2.

FDTD Method

The FDTD method developed by Yee²⁴ is a direct numerical solution of Maxwell’s curl equations in the time domain. Maxwell’s curl equations are discretized in space and time, resulting in six update equations for each Cartesian component of the electric and magnetic fields. These equations are updated in a leapfrogging manner as the simulation steps through time. The problem space is discretized into cubic voxels, each about $\lambda ∕20$ in size, to ensure the stability of the computational results. As a result, the computational complexity of a 3-D problem rapidly increases with problem size relative to wavelength. Satisfying memory requirements to simulate a scattering problem of realistic size quickly becomes impossible using typical desktop resources. However, due to the regular structure of the simulation grid, the problem is easily parallelized across many compute nodes to make large problems tractable when using a parallel computing cluster.

The same FDTD code used for a previous study²² was used for the superposition study. It uses the scattered-field-only formulation^{25, 26} to minimize memory requirements for storage of electric and magnetic field components. Due to this formulation, the incident $x$ -polarized plane wave is computed analytically and does not propagate via the FDTD update equations. The incident field is therefore immune from phase errors that may occur due to numerical dispersion often associated²⁶ with FDTD. The FDTD grid is terminated using Berenger’s split-field perfectly matched layer (PML) absorbing boundary condition.²⁷ The PML condition has been used extensively in FDTD simulations of optical scattering.

3.

Superposition Technique

The FDTD method computes only the electromagnetic field values in the near field of the problem space. To obtain the far-field scattering pattern, a near-to-far-field transform must be employed.^{26, 28} Using the discrete Fourier transform with the origin being the center of the problem space, the time-domain field values from the FDTD simulations over a surface $S$ were converted to frequency-domain surface current densities. The surface $S$ was chosen to be a cube, 20 grid voxels smaller than the total problem space size on all sides. The entire scattering geometry was still contained within this reduced volume. The far-field electric and magnetic vector potentials $\mathbf{A}$ and $\mathbf{F}$ are calculated from these surface current densities and then used to determine the electromagnetic far-field values. The far-field scattering pattern is then determined from these electric and magnetic field values.

The superposition method presented here is a modified version of the near-to-far-field transform. Smaller simulations that require fewer computational resources may be linearly superimposed in the same transform to estimate the far-field scattering from a larger scattering problem. In the case of optical scattering by tissue, the subgeometry can be chosen to be a single cell, but multicell subproblems can also be used. Scattering information from several subgeometries can be superimposed in a manner so as to approximate the scattering from many cells arranged in a tissue. The generalized formulation for the far-field vector potentials is

The accuracy of the superposition estimate is determined by comparing the far-field scattering pattern predicted by the superposition method with the scattered far-field of the large problem. Because of the large variations with respect to scattering angle in the amplitude of the scattered field signal, a correlation factor $\left(\rho \right)$ of the superposition pattern is computed with respect to the true scattering pattern from the large problem geometry²⁹:

Because a perfect replication of the true signal would result in $\rho =1$ , the error of the superposition far field estimate can be defined as

The accuracy of the superposition algorithm in estimating the true far-field scattering from a larger scattering geometry such as a tissue is dependent on the degree of multiple scattering between each of the component geometries, such as single cells. This is because the superposition estimation assumes no mutual scattering interaction between the individual subgeometries being superimposed on one another. In the problems investigated, the incident plane wave travels in the $+z$ direction. Because biological cells are individually very weak scatterers at optical wavelengths, arranging cells parallel to the $z$ axis can result in significant scattering interaction. This will make a superposition case using single cells less accurate in approximating the far-field scattering from such an arrangement. A multicell structure arranged perpendicular to the direction of plane wave propagation, however, should have very small multiple scattering occur between adjacent cells, making a superposition approximation using single cells as the component geometry much more accurate relative to the parallel configuration.

For larger multicell scattering problems with multiple cells along each Cartesian axis, it is necessary to determine the minimum number of cells that must be included in each subproblem to achieve a desired level of accuracy in far-field estimation. This will depend greatly on the level of multiple scattering occurring between various parts of the geometry. At a minimum, the subgeometries should likely consist of groups of scatterers arranged along the axis of the plane wave propagation to include some information about scattering in that direction.

4.

Cell Construction Parameters

A cell construction script used in another paper²² was adapted for the creation of problem spaces containing multiple randomized complex cellular scatterers. For more accurate modeling of scattering by cells embedded in tissue, the background is matched to that of cytoplasm [ $n=1.36$ (Ref. 30)]. The refractive index values used for the nuclei and mitochondria were $n=1.40$ (Ref. 30) and $n=1.38$ (Ref. 31), respectively. The cuboidal cells shown in Fig. 1 have major diameters of $15\mu \mathrm{m}$ and minor diameters of $13\mu \mathrm{m}$ . However, due to computational constraints, the simulation was unable to accurately resolve a cell membrane, so no boundary exists at these diameters. Because other investigations using FDTD to investigate cell scattering have found that the membrane has a relatively small contribution to far-field scattering,¹⁵ this was considered a reasonable decision. The major and minor diameters of the heterogeneous ellipsoidal nuclei are 6 and $5\mu \mathrm{m}$ , respectively. The overlapping spatial inhomogeneities within each nucleus are $1\mu \mathrm{m}$ in diameter, and have randomly assigned refractive index values uniformly distributed in the range of $1.40\pm 0.03$ . The refractive index in the overlapping regions is assigned in sequential order of sphere placement. The number of spherical inhomogeneities was chosen to equal roughly 10 times the volume of each nucleus to ensure heterogeneity throughout each nucleus. The organelles within each cell are divided evenly by volume into two groups: one ellipsoidal group with major and minor diameters of 1.5 and $0.5\mu \mathrm{m}$ and another spherical group with diameters of $0.5\mu \mathrm{m}$ . The location and orientation of cellular components are chosen via random number generation and organelle placement is constrained so they are nonoverlapping (see Fig. 2 ). The numeric seed for the random number generator can be fixed to enable a particular geometry to be repeated, or it may be set to change to enable original cell construction and to test for sensitivity of results to particular cellular configurations.

Fig. 1

Rendering of eight-cell cube with randomly placed homogeneous nuclei and organelles.

Fig. 2

Refractive index map of a 2-D slice of a 3-D 27-cell problem geometry with heterogeneous nuclei. The overlapping spherical heterogeneities are seen within the nuclei and the organelles within the cell are constrained to be nonoverlapping.

5.

Results

5.1.

Eight Cells

Eight different randomly generated cells with heterogeneous nuclei and 6% organelle density by volume were arranged into a $2\times 2\times 2$ cube (see Fig. 1 for the arrangement of cells). An FDTD simulation was performed to calculate the scattering of a $+z$ traveling plane wave by this geometry. A near-to-far-field transform²⁸ was used to obtain the far-field scattering from a simulation of all eight cells. The geometry was then split into two subgroups, each having four cells arranged in a $2\times 2$ square in the $yz$ plane (or $1\times 2\times 2$ geometry). A similar subgrouping was made using two $2\times 2$ squares each arranged in the $xy$ plane (or $2\times 2\times 1$ geometry). Smaller subdivisions of the geometry were also made by grouping the cells in four subgroups of size $2\times 1$ arranged in either the $y$ or $z$ direction, and finally each of the eight cells was treated as a subgroup on its own. Each of these subgroups was simulated separately using FDTD and near-field scattering data were recorded. The scattering data from each subgroup case were then assembled by the superposition algorithm into a scattering approximation of the original $2\times 2\times 2$ cube of cells.

Figure 3 shows the far-field scattering pattern $F_s\left(\theta \right)$ averaged over all $\varphi$ angles for the original eight cell geometry along with several results of the superposition method using various different subgeometries. While the scattering was calculated for all $\theta$ , only $0⩽\theta ⩽30$ is shown here for easier differentiation between traces. The legend entry $2\times 2\times 2$ indicates the original eight cell geometry results, while, for example, the entry $1\times 1\times 2$ indicates the results of the superposition of the far-field scattering of four separate groups of two cells each arranged along the $z$ axis. Additionally, Fig. 3 shows the calculated scattering cross section and anisotropy resulting from the far-field pattern in each case.

Fig. 3

Average scattering $F_s\left(\theta \right)$ over all $\varphi$ angles along with scattering cross section and anisotropy of the $2\times 2\times 2$ different cell geometry and superposition results for various subgeometries.

The superpositions using subgroups $1\times 2\times 2$ and $1\times 1\times 2$ show the greatest agreement with the scattering from the original geometry, with errors of 0.005 and 0.145, respectively. The remaining superposition estimates have error figures that are almost exactly 1, which indicates that there is barely any correlation between those superposition estimates and the true scattered field. This should not be surprising as the incident plane wave is traveling in the $+z$ direction and the cellular components are all relatively weak scatterers compared to the cytoplasm background. This will cause most multiple scattering between cells to occur roughly along the $z$ axis. Those subgroups of cells whose configurations are along the $z$ axis will preserve most effects of this multiple scattering, while those that are arranged only along the $x$ the or $y$ axis will effectively ignore them.

The results shown in Fig. 4 demonstrate the effectiveness of the superposition method in predicting scattering for multiple cells arranged in a plane perpendicular to the direction of incident light. Because the cellular geometry is $2\times 2\times 1$ , there is no multiple scattering between cells along the $z$ axis. The superposition of $1\times 2\times 1$ geometries is actually slightly less accurate with an error of 0.059 compared with an error of 0.002 of the superposition of single cells $(1\times 1\times 1)$ , their differences in error are relatively small. Because there is much less scattering interaction between the cells in this case, it can be seen as an example of the type of problem that would greatly benefit from the superposition of many small scattering simulations, such as single cell simulations.

Fig. 4

Average scattering $F_s\left(\theta \right)$ over all $\varphi$ angles along with scattering cross section and anisotropy of the $2\times 2\times 1$ different cell geometry and superposition results for various subgeometries.

While the problem geometry and superposition results shown in Figs. 1 and 3 show the eight cell problem using eight differently constructed cells, similar results were seen in the case where the exact same cell is replicated eight times into a $2\times 2\times 2$ cube. This demonstrates that the superposition technique can work equally well for replicating the same cell into a larger geometry or assembling different subgroups into a larger arrangement.

5.2.

Results from Increasingly Complex Scattering Geometries

Expanding the size of the previous scattering geometry in the two directions perpendicular to the light propagation axis, a $3\times 3\times 2$ cell geometry was created [see Fig. 5 ]. The cells in this case are generated in the same manner as those previously seen. The far-field and superposition scattering pattern results for this geometry are seen in Fig. 6 . Little qualitative changes are seen in the results when compared to the $2\times 2\times 2$ cell geometry, with the $z$ arranged subgeometries resulting in the most accurate far-field scattering predictions. The errors for the $1\times 3\times 2$ and the $1\times 1\times 2$ superposition cases are 0.104 and 0.231, respectively.

Fig. 5

(a) Simplified diagram of $3\times 3\times 2$ cell scattering geometry and (b) one possible subdivision into $3\times 3\times 1$ subproblems.

Fig. 6

Average scattering over all $\varphi$ angles along with scattering cross section and anisotropy of the $3\times 3\times 2$ different cell geometry and superposition results for various subgeometries.

When a third layer of cells is added in the $z$ dimension [see Fig. 7 ], however, some of the trends seen in previous superposition results change. Each cell in this $3\times 3\times 3$ case is again generated with a different random seed and each has a heterogeneous nucleus. The importance of using $z$ arranged subgeometries is still observed in the far-field scattering pattern of this larger problem, as seen in Fig. 8 . The patterns show that the $1\times 3\times 3$ and $1\times 1\times 3$ superposition cases remain in the best agreement with the scattering pattern of the large problem, with errors of 0.011 and 0.227, respectively. Interestingly, the errors of the other superposition cases are no longer nearly 1, with all values of the remaining cases being around 0.88. While these estimates still have quite low correlation with the true far-field pattern, the correlation is no longer close to zero.

Fig. 7

(a) Simplified diagram of $3\times 3\times 3$ cell scattering geometry and (b) one possible subdivision into $3\times 3\times 1$ subproblems.

Fig. 8

Average scattering $F_s\left(\theta \right)$ over all $\varphi$ angles along with scattering cross section and anisotropy of the $3\times 3\times 3$ different cell geometry and superposition results for various subgeometries.

The anisotropy and scattering cross section estimates from any of the multicell superposition cases, however, are in much closer agreement with those of the $3\times 3\times 3$ problem than was true in the $2\times 2\times 2$ or $3\times 3\times 2$ cell cases. This appears to be a result of the choice of cell construction parameters and cell arrangement in the $3\times 3\times 3$ problem geometry. Multiple scattering from the front layers of cells actually has a shading effect on the rear layers, making those rear layers have a decreasing impact on the total amount of light scattered by the volume.³² In certain subdivision cases, however, the superposition rule assumes that all the superimposed layers are exposed to the same incident light intensity, causing their scattering cross sections to add nearly linearly (discounting destructive interference due to phase differences). In the case where three $3\times 3\times 1$ subgroups are superimposed along the $z$ axis, for example [see Figs. 5 and 7], the resulting prediction for scattering cross section is roughly a sum of the cross sections of the subgroups. It happens that the cross section of each of these subgroups is nearly one third of that of the full $3\times 3\times 3$ problem. Similar phenomena are seen for smaller subgroups that are not arranged along the $z$ axis. Thus, in this one specific case the superposition approximation provides an accurate estimator for cross section and even anisotropy regardless of the way the problem is subdivided. This level of accuracy when using any choice of subgroup is not necessarily repeatable if different cells are used in the same arrangement, or even if the same cells are used but spaced slightly differently. Thus, the agreement between, for example, the $3\times 3\times 1$ subdivision case prediction for scattering cross section and the $3\times 3\times 3$ cross section does not represent a trend toward any superposition case being a good predictor of scattering cross section. In addition, the far-field scattering pattern prediction of the various superposition cases follows the same trend as previous cases, with the $1\times 3\times 3$ and $1\times 1\times 3$ predictions being much more accurate than any other choice of subgeometry. However, accurate superposition results are still obtained in these cases when a suitable $z$ arranged subgeometry is used.

To determine if the observed trends in scattering cross section prediction and error measurements continued with increasing tissue thickness, an even larger simulation consisting of 36 cells in a $3\times 3\times 4$ configuration was simulated. The results are shown in Fig. 9 . Once again, the subgroups arranged along the $z$ axis resulted in the best agreement with the large simulation results. The error figures in the $1\times 3\times 4$ and $1\times 1\times 4$ cases were 0.044 and 0.239, respectively, while the remaining error figures were all unity. This suggests that the optimal method for choosing the superposition subgrouping will remain the same for larger problems. Looking at the bar graphs in Fig. 9, one notices that the $3\times 3\times 1$ and $1\times 3\times 1$ superpositions have actually begun to overpredict the scattering cross section of the true problem.

Fig. 9

Average scattering $F_s\left(\theta \right)$ over all $\varphi$ angles along with scattering cross section and anisotropy of the $3\times 3\times 4$ different cell geometry and superposition results for various subgeometries.

For this particular choice of cell geometry and arrangement, this $3\times 3\times 3$ configuration appears to be a crossover point in scattering cross section prediction when using the $3\times 3\times 1$ subdivision as the number of cells along the $z$ direction of the scattering geometry are increased. In the case where the scattering geometry is only two cells thick in the $z$ direction, poor choices in superimposed subgroups, such as the one seen in Fig. 5, cause an underprediction in scattering cross section. Because of the diminishing real contribution to the overall scattering cross section by the third layer of cells, that same poor superposition choice, seen in Fig. 7, is much more likely to result in an overprediction of the scattering cross section in the $3\times 3\times 3$ geometry. Adding an additional layer of scattering cells to the geometry along the $z$ direction causes these same superposition arrangements to further overpredict the scattering cross section of the large problem. Nonetheless, the scattering pattern plots in Figs. 8 and 9, as well as their corresponding error figures $\left(\rho \right)$ clearly demonstrate that the rules on choosing an optimal subgeometry learned from previous problems remain in effect, even if one is only interested in the figures of merit such as scattering cross section or anisotropy rather that the far-field scattering pattern. Multiple scattering along the propagation axis remains the most important factor in the choice of superposition geometries.

6.

Discussion

Scattering from single cells or a small number of cells has been extensively investigated^{9, 15, 18, 20} computationally using FDTD. While these studies remain useful for determining the effects of cell morphology on light scattering and propagation, the scattering from larger tissue structures containing many cells may be of more practical interest to researchers whose instrumentation measures light scattered from relatively large tissue sections. The results of multicell simulations such as those shown here suggest that using computational methods for bulk scattering parameter prediction requires multiple cells to be included in the problem geometry.²² In addition, the consistency of the far-field scattering pattern between the $3\times 3\times 2$ , $3\times 3\times 3$ , and $3\times 3\times 4$ geometries suggest that the number of cells required for an accurate prediction may not be much larger than is seen here. The number is also likely to be tissue-dependent. Regardless of the exact number of cells required to obtain an accurate bulk parameter prediction, that number is likely to be sufficiently large to require significant computational resources.

Computational requirements for FDTD simulations of such large problems render them intractable using most current resources. Using the near-field scattering results of a single row or a small number of rows of cells aligned with the direction of plane wave propagation, the linear superposition method already outlined could be used to provide an estimate of far-field scattering by a large tissue using practical computational resources. Using this approach to estimate a larger scattering problem effectively reduces the dimensionality of the problem (for example, $3\times 3\times 3$ cells to $1\times 1\times 3$ cells), reducing the processing time and drastically reducing the memory requirements to perform the FDTD simulation. In the case of $3\times 3\times 3$ differently generated cells estimated by nine separate $1\times 1\times 3$ cell simulations, the computation time would be reduced from 5300 to 4500 CPU hours and the maximum memory requirement would be $12\mathrm{Gbytes}$ versus $130\mathrm{Gbytes}$ for the full $3\times 3\times 3$ problem.

Some care must be exercised as to how a large scattering geometry is broken up to maintain a relatively low error in far-field scattering predictions by the superposition algorithm. At a minimum, one should choose a subgeometry that extends along the full length of the scattering problem in parallel with the light propagation axis. This accounts for multiple scattering between cells along the propagation axis, which are the bulk of the light scattering effects in a small tissue. This choice of subgeometries is valid only for mostly forward-scattering problems such as biological tissue. For more densely scattering geometries, or those with higher refractive indices relative to the background, it may be necessary to simulate a subsection of the geometry and examine the resulting far-field scattering pattern. As an example, a $z$ directed column, much like a $1\times 1\times 3$ cell case, could be used as a test geometry. If the majority of the scattered light is contained near $\theta =0$ , the geometry is a likely candidate for the superposition method. The exact choice of an acceptable angular region is dependent on the depth in $z$ of the large problem being approximated by the superposition method.

Selecting a suitable subgeometry to preserve the majority of multiple scattering interaction along the $z$ axis is obviously much simpler when the cells are arranged into a rectangular grid, as in the preceding geometries. However, this kind of cellular arrangement may not be ideal for describing certain tissue types. It was also found that this regular arrangement of cells was also quite sensitive to small changes in $z$ separation between successive cell layers. While these changes did not greatly affect the majority of the far-field scattering pattern, small changes in separation distances between cells cause significant changes to the pattern at very small angles, greatly affecting the scattering cross section and anisotropy results. The corresponding superposition results using subgeometries that extend along the $z$ axis, however, continued to accurately predict this behavior.

While this method is demonstrated here as a way to reduce the problem size of a large tissue-like scattering problem, it can also be seen as a method to use single-cell scattering data to simulate a periodic structure of individual cells, such as the single-file queue in a flow cytometer.¹⁸ The superposition method could be used as an alternative to a periodic boundary condition in the lateral direction of the queue, mitigating the need for a custom FDTD boundary setting for this kind of problem.

The use of other cellular arrangements does not necessarily prove to be an impediment to using the superposition method. In the preceding problems, a logical choice was to divide the geometry along the cell boundaries, but there is nothing preventing other subdivision methods. So long as the divisions do not go through individual scatterers, such as nuclei or organelles, any subdivision that extends the length of the problem along the propagation axis could work. This is provided that this subdivision is still wide enough in the $x$ and $y$ dimensions to allow for most of the energy that is scattered at small angles to propagate to the rear of the subgeometry. This will maintain the accuracy of the approximation of no scattering interaction between subgeometries.

If the surface current densities resulting from the FDTD simulations, ${\mathbf{J}}_s$ and ${\mathbf{M}}_s$ , are saved in a data file, they may be reused by any subsequent superposition calculation to piece together any number of far-field scattering predictions without incurring any computational costs save those for the discrete Fourier transforms in the superposition method. This can allow enable scientists to compile a database of FDTD surface current data from any number of scattering geometries and use the superposition algorithm to make predictions on scattering by any spatial combination of those geometries. The spatial configuration of the geometries will, of course, determine the accuracy of the superposition prediction, with arrangement of the scatterers in a plane perpendicular to the direction of plane wave propagation being optimal.

7.

Conclusion

A coherent linear superposition method for estimating the plane wave far-field scattering pattern from multiple biological cells computed by the FDTD method was presented. Using the FDTD-simulated near-field scattering results from small numbers of complex scatterers, an estimate of the far-field pattern from a large group of those same scatterers can be made. This method can be used to reduce the computational cost of some FDTD simulations by enabling a single large scattering problem to be broken into smaller FDTD problems with more practical computational requirements. The results of some example cases demonstrated certain considerations that must be made for choosing the minimum number and arrangement of cells to simulate in large multicell scattering problems. In particular, it was found that the method works best in cases where there is very little multiple interaction between adjacent scatterering subgroups, so the far-field pattern can be estimated as a superposition of the scattering from small numbers of cells.