1. Introduction

The interaction between electromagnetic waves and tissue is important for the development of diagnostic and therapeutic applications of light in medicine. For example, scattering of light within tissue simultaneously makes possible and limits high-resolution imaging techniques such as reflectance-mode confocal microscopy [1] and optical coherence tomography (OCT) [2]. To better understand the interaction of light within tissue and, thus, the possibilities of these and other imaging modalities, it is useful to have a model of tissue from which the optical properties can be derived.

Currently, several models exist with which to simulate the optical properties of tissue. The finite-difference time-domain (FDTD) method has been employed to simulate scattering from single cells by direct solution of Maxwell’s equations [3,4]. For bulk tissue, Mie theory [5,6] has been utilized by approximating the tissue with a fractal distribution of spheres [7, 8]. A fractal-continuous distribution [9, 10] has been used to investigate light scattering in tissue based on the sample refractive-index correlation function.

With the use of phase-contrast microscopy, it has been shown that the scattering centers in biological tissue aggregate into complicated forms similar to those of fractal objects [11, 12]. A model based on a fractal size distribution of spherical scatterers is attractive in that it provides an approximation to the complicated spatial refractive-index variations within bulk tissue. Therefore, it provides the ability to simulate beam propagation through tissue and, thus, to describe models of imaging modalities that rely on its scattering properties [13]. The optical properties determined from the tissue model are also useful for Monte Carlo simulations of light transport in tissue [14].

The discrete particle model of biological tissue [8] is an attractive and physically appealing approach, however, it lacks a rigorous treatment of multiple and correlated scattering effects. Multiple scattering refers to the situation in which the scattered wave from each particle contributes to the wave incident upon the other scatterers [15], whereas correlated scattering refers to the interference effects between the scattered fields from all particles [16]. In this paper, we define multiparticle scattering to include both multiple and correlated scattering effects. The model’s use of a wavelength-independent packing factor [17, 18] is sufficient for describing correlated scattering among the small-sized particles of the fractal distribution that scatter in the Rayleigh regime. However, the distribution contains important contributions from Mie scatterers and their density makes the likelihood of correlated scattering, and especially multiple scattering, among these particles high. In addition, the effect of multiparticle scattering upon the anisotropy of the discrete particle model has not been investigated to date. Using a field-based computational method that extends Mie theory to the multisphere case, we simulated the optical properties of a fractal aggregate, further developing the discrete particle model by incorporating the effects of multiparticle light scattering. We show that multiparticle scattering has important effects upon the optical parameters simulated, which include the scattering coefficient, anisotropy, reduced scattering coefficient, and the phase function. In addition, we demonstrate the importance of incorporating an enhanced contribution of Rayleigh scattering to the model for reasonable agreement with measured data from skin tissue (dermis).

2. Theory

Following Schmitt and Kumar [8], we describe the scatterers within biological tissue in terms of a distribution of spherical particles, of varying sizes but equal, uniform refractive indices within a background medium (of constant refractive index). The scatterer size distribution for dense fibrous tissue, such as the dermal layer of skin tissue, follows a fractal power law distribution [12]. If the sphere diameter is denoted d, the volume fraction density function in terms of the sphere size is defined as:

so that the fraction of the total volume occupied by spheres in the range [d, d + Δd] is equal to η(d)Δd, D_f is the fractal dimension, K ₁ is a scaling constant, d _min is the smallest sphere size, and d _max is the largest sphere size of the distribution. The total volume fraction of the aggregate F _v is defined by the integral F _v = ∫^∞ ₀(d) dd.

The variations in the refractive index of the model are determined using weighted averages of typical refractive-index values for the constituents of biological tissue, to obtain background and particle refractive indices, n¯_bkg and n¯_p, respectively. For soft biological tissue, the background refractive index has been estimated to be approximately 1.35 [8], calculated by estimating the volume fractions of intracellular and interstitial fluids present in tissue. The particle refractive index is determined by calculating a quantity known at the mean index variation (Δn), based on the refractive indices and relative volume fractions of the collagen fiber, nuclei and organelle constituents of dermis. This quantity has been estimated to be equal to 0.10 [8], and the particle refractive index is given by n¯_p= n¯_bkg + (Δn) = 1.45.

To be able to simulate the optical properties of the model, an aggregate was created by randomly distributing spheres throughout a spherical region of volume V. (A sphere is defined to be part of the aggregate if its center is located within this region.) The volume V was chosen to be suitably large so that the optical properties computed by Mie theory sufficiently converge to those of the infinitely large aggregate. This also ensures that the statistical properties of bulk tissue are approximated as well as possible, given the computational constraints of the simulation. The spheres corresponding to the smallest particle size of the distribution, that scatter light within the Rayleigh regime, were excluded from the aggregate due to the computational resource limitations of the simulation, however, their contribution to the total optical properties was computed separately, as shown later. Similarly to Schmitt and Kumar [8], we chose a total of 10 discrete sphere sizes. However, in contrast to their choice of diameters equally spaced on a logarithmic scale, we sample sphere sizes linearly to best retain the fractal properties of the distribution. We also chose particle sizes ranging from 50 nm to 20 μm; the contribution of the particle sizes outside this range to the total optical properties of the model is negligible. The discrete particle volume fractions were weighted according to Eq. 1, so that their sum, the total volume fraction F _v, was equal to 0.2. Examples of the aggregates created are shown in Fig. 1 and show the higher density of larger particles present in the aggregate with a lower fractal dimension. The sizes of the two aggregates are different, since they have different particle densities but the number of spheres in each are approximately equal. For ease of comparison, they are displayed in Fig. 1 with the same axis scale.

We employ a solution to the multisphere scattering problem provided by Xu [19-23], known as the Generalized Multiparticle Mie solution (GMM) [24]. To model interaction (and, thus, multiple scattering) effects between spheres, this method expands the scattered fields in terms of vector spherical wave functions (VSWFs) that impinge on all other spheres within the aggregate. By also expanding the incident wave in terms of VSWFs, and making use of addition theorems for these functions, multiple scattering effects can be solved analytically. By taking account of the phase relation and direction of all scattered waves, correlated scattering in the far-field is solved by computing their interference. Inputs into the GMM FORTRAN program [23] are the spatial locations, size, and relative refractive indices of all spheres, as well as the wavelength of the unpolarized incident plane wave. The output optical parameters include the total scattering cross-section σ_s,agg, anisotropy g _agg, and the phase function P _agg averaged over all azimuthal angles. The total scattering coefficient of the aggregate may then be calculated as:

 

Fig. 1. Visual representation of example aggregates (excluding the smallest sphere size of 50 nm) for a fractal dimension of 3.5 and 3.7 shown in (a) and (b), respectively (scale in micrometers).

Download Full Size | PDF

For weakly-scattering subwavelength-diameter spheres, multiparticle scattering is dominated by correlated scattering which may be described by the use of a packing factor [17, 18]. This reduces the bulk optical cross-section of these scatterers so that their effective volume fraction density function can be written as:

By using the correction factor implicit in this equation, the contribution to the total scattering cross-section from the smallest particle size initially excluded from the aggregate can be computed and included a posteriori. Scattering from skin tissue typically includes a strong Rayleigh contribution [25-27], much greater than that predicted by the fractal distribution alone. As seen in Eq. 3, the use of this correction factor does not change the wavelength dependence of the Rayleigh scattering and, therefore, an added Rayleigh scattering component to the scattering coefficient can be written without explicitly accounting for the scatterer sizes as:

where K ₂ is a parameter that varies the magnitude of this additional component. We will assume that correlated scattering among Rayleigh scatterers does not significantly change the overall isotropic scattering behavior exhibited by these particles and, thus, the effect of their contribution upon the total anisotropy can be derived as:

The reduced scattering coefficient, defined as μ′_s = (1 - g)μ_s, is therefore given by the expression:

The phase function is derived from the angular scattered irradiance (per unit incident irradiance) of a medium. It is a probability distribution function for the scattering angle of a photon relative to its incident direction. If multiparticle scattering effects are ignored, the volume-averaged phase function for M particle sizes is given by:

which is normalized so that ∫^π ₀ P(θ)2π sin θ dθ = 1. The phase function for Rayleigh scatterers (given unpolarized incident light) may be written as $P_{\mathrm{Ray}}\left(\theta \right)=\frac{3}{16\pi }\left(1+{\cos }^2\theta \right),$, so that the total phase function of the aggregate can be shown to be equal to:

The wavelength range simulated was between 400 and 1800 nm, incorporating both the “therapeutic window” and part of the near-infrared region of the electromagnetic spectrum. The wavelength increment of the simulation was chosen to be 200 nm, due to the substantial computation time required by the program. (For a wavelength of 400 nm, a CPU time of approximately 160 hours was necessary to simulate an aggregate.) Larger incident wavelengths require less CPU time due to the lower order of the incident and scattered field expansions necessary for convergence. A total of five random realizations of the aggregate corresponding to a particular fractal dimension were created and the optical properties averaged for each data point. In computing the optical properties without taking into account multiparticle scattering effects, we chose a wavelength increment of 1 nm due to the substantially smaller computation time necessary.

3. Results and discussion

The wavelength dependence of the scattering coefficient, anisotropy, and reduced scattering coefficient from the simulations of several fractal aggregates is shown in Fig. 2. The standard error of the mean for each data point was not more than 0.35 mm ^-1 for the scattering coefficient and 0.002 for the anisotropy. Based on spatial frequency analysis, the fractal dimension for a particular soft tissue sample has been estimated to be around 3.6 [11,12]. We chose three fractal dimensions in the range 3.5 ≤ D_f ≥ 3.7, which lead to optical properties consistent with measured data for skin tissue.

The wavelength dependence of the scattering coefficient, not taking into account multiparticle scattering effects, can be approximated by the power-law relation μ_s ∝ = λ^3-D_f. This differs from the relation reported by Schmitt and Kumar [8], which has an exponent of 2 - D_f. The discrepancy is the result of our sampling of the discrete particle sizes on a linear scale (as opposed to a logarithmic scale) to represent a fractal distribution. Our result is consistent with the relation reported by Wang [28], in which the sphere sizes were also sampled linearly. This simple relation does not hold, however, if multiparticle scattering is taken into account, and for all fractal dimensions, it can be seen that the apparent exponent of the scattering coefficient wavelength dependence is reduced. The reduction in magnitude of the scattering coefficient is most significant at lower wavelengths, at which multiparticle scattering effects associated with larger Mie scatterers dominate.

 

Fig. 2. Anisotropy (a) and scattering coefficient (b) for multiparticle (solid lines) and non-multiparticle (dashed lines) scattering effects. The correspondence between fractal dimension and line color is indicated in the legend. The dotted lines in (a) shows the effect of the Rayleigh contribution upon the anisotropy. The reduced scattering coefficient is shown in (c) for multiparticle (solid lines) and an added contribution of Rayleigh (dotted lines) scattering. (No Rayleigh contribution is included for the aggregate of fractal dimension 3.7.) Empirical data points (shown as markers) are from [27] and the references therein [32-36].

Download Full Size | PDF

In general, the anisotropy of the aggregates is decreased due to multiparticle scattering; for wavelengths above around 1000 nm, the effect is greatest for aggregates with lower fractal dimensions. For wavelengths below this value, the contribution from small Rayleigh scatterers becomes prominent due to the λ_-4 dependence in Eq. 4, which for aggregates with larger fractal dimensions, and hence larger relative volume fractions of small scatterers, reduces the anisotropy. This strong reduction in anisotropy for shorter wavelengths, which is not present in the discrete particle model of Ref. [8], is consistent with the empirical fit described by van Gemert et al. [29] based on the measured data of Bruls and van der Leun [30], and that reported by Beek et al. [31].

The reduction in the scattering coefficient is accompanied by a corresponding reduction in the anisotropy; these effects partially cancel when calculating the reduced scattering coefficient. Inspection of the experimentally measured data sets in [27, 32, 33] shows that the exponent of the reduced scattering coefficient for shorter wavelengths approaches the Rayleigh limit of λ^-4, which is not accounted for in the fractal model since there is typically an underestimation of the number of particles that scatter in the Rayleigh regime. (Multiparticle scattering within the aggregate reduces this exponent even further.) This demonstrates that, for human dermis, an added Rayleigh scattering contribution to the discrete particle model will generally substantially improve its agreement with empirical data at low wavelengths. The value of K ₂, which sets the magnitude of this contribution, was chosen by inspection so that the model gives reasonable results when compared to measured data for both the reduced scattering coefficient and anisotropy. For the wavelength defined in micrometers, the value of K ₂ was set to 0.07 and 0.03 for the aggregates with a fractal dimension of 3.5 and 3.6, respectively. The result of this addition of Rayleigh scattering to these aggregates is shown in Fig. 2. There is reasonable agreement between the reduced scattering coefficient of the model and that of experimental data for human and porcine dermis, both of which are similar in structure [34, 37].

 

Fig. 3. Phase functions of the aggregate for a wavelength of 633 nm and 1300 nm shown in (a) and (b), respectively, with a fractal dimension of 3.6. The dashed lines (blue) and solid lines (red) represent non-multiparticle and multiparticle scattering, respectively. The dotted line (green) represents multiparticle scattering including a contribution of Rayleigh scattering. The inset shows the phase function in the range from 0 to 15 degrees.

Download Full Size | PDF

The shape and accuracy of the phase function can significantly impact the results of Monte Carlo simulations of light transport in turbid media [38, 39]. The phase function for a specific realization of the fractal aggregate is shown in Fig. 3 for wavelengths of 633 and 1300 nm. For both wavelengths, there is a significant increase in the probability of direct forward scattering (corresponding to a scattering angle θ = 0°), due to the fact that the majority of the scattered waves are in phase for this particular scattering direction [40,41] and, hence, is present for every realization of the aggregate. Even though this forward peak increases the anisotropy slightly, multiple scattering effects result in an overall reduction in anisotropy. The specific locations of the local maxima and minima that form the high-frequency structure in the phase function are dependent on the random positions of the spheres within the aggregate. These strong resonances, seen in Fig. 3, are not likely to be present in real tissues due to the averaging effects of the non-spherical scatterers and random, continuous variations in refractive index. Averaging the phase function of many realizations of the aggregate, however, reduces the fluctuations caused by these resonances and, thus, we expect our results to be applicable to simulating the effects of real tissue. Additional Rayleigh scatterers increase the probability of direct backscattering (corresponding to θ = 180°), which, as expected, is higher for lower wavelengths.

The close-packed nature of scatterers in soft biological tissue can sometimes be overlooked, particularly when high-index contrast microspheres are used in experiments as tissue phantoms [42]. To achieve realistic optical properties, such phantoms have very small volume fractions and, thus, we expect multiparticle scattering effects to be negligible in such phantoms. The effects of multiparticle scattering are, broadly, to reduce the magnitude and exponent of the wavelength dependence of the scattering coefficient, and to reduce the anisotropy. Comparison with experimental data, however, demonstrates the importance of also adequately accounting for an additional Rayleigh scattering contribution.

4. Conclusion

We have shown that, when taken into account, multiparticle scattering has a significant effect upon the optical properties of a fractal soft tissue model. Combined with a Rayleigh scattering contribution, the wavelength dependency of the resulting reduced scattering coefficient is in reasonable agreement with experimentally determined values for skin tissue (dermis). A forward peak in the phase function due to correlated scattering is observed, however the effect of multiple scattering reduces the total anisotropy over all wavelengths simulated. The model presented here, or variations of it to account for different tissue types, will provide a more accurate basis for a range of tissue optics problems. Further to the results presented here, it is possible to determine the Muüller matrix, which relates the incident and scattered Stokes parameters, from the amplitude-scattering matrix computed from the simulation. Thus, it would be possible to model the depolarization of light propagating through the skin tissue model.