We rigorously account for the effects of multiparticle light scattering from a fractal sphere aggregate in order to simulate the optical properties of a soft biological tissue, human skin. Using a computational method that extends Mie theory to the multisphere case, we show that multiparticle scattering significantly affects the computed optical properties, resulting in a reduction in both scattering coefficient and anisotropy for the wavelengths simulated, as well as a significantly enhanced forward peak in the simulated phase function. The model is extended to incorporate the contribution of Rayleigh scatterers, which we show is required to obtain reasonable agreement with experimentally measured optical properties of skin tissue.
© 2007 Optical Society of America
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  and optical coherence tomography (OCT) . 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 . The optical properties determined from the tissue model are also useful for Monte Carlo simulations of light transport in tissue .
The discrete particle model of biological tissue  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 , whereas correlated scattering refers to the interference effects between the scattered fields from all particles . 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).
Following Schmitt and Kumar , 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 . 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, Df is the fractal dimension, K 1 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 = ∫∞ 0(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 , 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 , 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 , 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) . 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  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:
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 2 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 ∫π 0 P(θ)2π sin θ dθ = 1. The phase function for Rayleigh scatterers (given unpolarized incident light) may be written as , 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 ≤ Df ≥ 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-Df. This differs from the relation reported by Schmitt and Kumar , which has an exponent of 2 - Df. 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 , 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.
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. , is consistent with the empirical fit described by van Gemert et al.  based on the measured data of Bruls and van der Leun , and that reported by Beek et al. .
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 2, 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 2 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].
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 . 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.
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.
DHPS is supported by the William and Marlene Schrader Trust. The computational simulations were run on an SGI Altix supercomputer provided by the ARRC computing facility at iVEC in Western Australia.
References and links
1. T. Wilson and C. J. R. Sheppard, Theory and practice of scanning optical microscopy (Academic Press, London,1984).
2. D. Huang, E. A. Swanson, C. P. Lin, J.S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254, 1178–1181 (1991). [CrossRef] [PubMed]
3. A. Dunn and R. Richards-Kortum, “Three-dimensional computation of light scattering from cells,” IEEE J. Sel. Topics Quantum Electron. 2, 898–905 (1996). [CrossRef]
4. R. Drezek, A. Dunn, and R. Richards-Kortum, “A pulsed finite-difference time-domain (FDTD) method for calculating light scattering from biological cells over broad wavelength ranges,” Opt. Express 6, 147–157 (2000), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-6-7-147. [CrossRef] [PubMed]
5. H. C. van de Hulst, Light scattering by small particles ( Dover Publications, New York,1981).
6. C. F. Bohren and D. R. Huffman, Absorption and scattering of light by small particles (Wiley-Interscience, New York,1983).
7. B. Gélébart, E. Tinet, J. M. Tualle, and S. Avrillier, “Phase function simulation in tissue phantoms: a fractal approach,” Pure Appl. Opt. 5, 377–388 (1996). [CrossRef]
8. J. M. Schmitt and G. Kumar, “Optical scattering properties of soft tissue: a discrete particle model,” Appl. Opt. 37, 2788–2797 (1998). [CrossRef]
10. C. J. R. Sheppard, “Fractal model of light scattering in biological tissue and cells,” Opt. Lett. 32, 142–144 (2007). [CrossRef]
12. G. Kumar and J. M. Schmitt, “Micro-optical properties of tissue,” in Advances in Laser and Light Spectroscopy to Diagnose Cancer and Other Diseases III: Optical Biopsy,R. R. Alfano, ed., Proc. SPIE 2679, 106–116 (1996). [CrossRef]
13. J. M. Schmitt and A. Knüttel, “Model of optical coherence tomography of heterogeneous tissue,” J. Opt. Soc. Am. A 14, 1231–1242 (1997). [CrossRef]
15. C. Smart, R. Jacobsen, M. Kerker, J. P. Kratohvil, and E. Matijević, “Experimental study of multiple light scattering,” J. Opt. Soc. Am. 55, 947–955 (1965). [CrossRef]
17. V. Twersky, “Acoustic bulk parameters in distributions of pair-correlated scatterers,” J. Acoust. Soc. Am. 64, 1710–1719 (1978). [CrossRef]
20. Y. -L. Xu, “Electromagnetic scattering by an aggregate of spheres: far field,” Appl. Opt. 36, 9496–9508 (1997). [CrossRef]
21. Y. -L. Xu, “Electromagnetic scattering by an aggregate of spheres: asymmetry parameter,” Phys. Lett. A 249, 30–36 (1998). [CrossRef]
22. Y. -L. Xu, “Electromagnetic scattering by an aggregate of spheres: theoretical and experimental study of the amplitude scattering matrix,” Phys. Rev. E 58, 3931–3948 (1998). [CrossRef]
23. Y. -L. Xu, gmm01s.f, http://atol.ucsd.edu/scatlib (2007).
24. Y. -L. Xu and B. Å. S. Gustafson, “A generalized multiparticle Mie-solution: further experimental verification,” J. Quant. Spectrosc. Radiat. Transfer 70, 395–419 (2001). [CrossRef]
26. R. Graaff, A. C. M. Dassel, M. H. Koelink, Mul F. F. M. de, J. G. Aarnoudse, and W. G. Zijlstra, “Optical properties of human dermis in vitro and in vivo,” Appl. Opt. 32, 435–447 (1993). [CrossRef] [PubMed]
27. A. N. Bashkatov, E.A. Genina, V. I. Kochubey, and V. V. Tuchin, “Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm,” J. Phys. D: Appl. Phys. 38, 2543–2555 (2005). [CrossRef]
28. R. K. Wang, “Modelling optical properties of soft tissue by fractal distribution of scatterers,” J. Mod. Opt. 47, 103–120 (2000).
31. J. F. Beek, P. Blokland, P. Posthumus, M. Aalders, J. W. Pickering, H. J. C. M. Sterenborg, and M. J. C. van Gemert, “In vitro double-integrating-sphere optical properties of tissues between 630 and 1064 nm,” Phys. Med. Biol. 42, 2255–2261 (1997). [CrossRef] [PubMed]
32. S. A. Prahl, “Light transport in tissue,” Ph.D. Dissertation, University of Texas at Austin (1988).
33. C. R. Simpson, M. Kohl, M. Essenpreis, and M. Cope, “Near-infrared optical properties of ex vivo human skin and subcutaneous tissues measured using the Monte Carlo inversion technique,” Phys. Med. Biol. 43, 2465–2478 (1998). [CrossRef] [PubMed]
36. E. K. Chan, B. Sorg, D. Protsenko, M. O’Neil, M. Motamedi, and A. J. Welch, “Effects of compression on soft tissue optical properties,” IEEE J. Quantum Electron. 2, 943–950 (1996). [CrossRef]
37. R. M. Lavker, G. Dong, P. Zheng, and G. F. Murphy, “Hairless micropig skin. A novel model for studies of cutaneous biology,” Am. J. Pathol. 138, 687–697 (1991). [PubMed]
38. J. R. Mourant, J. Boyer, A. H. Hielscher, and I. J. Bigio, “Influence of the scattering phase function on light transport measurements in turbid media performed with small source-detector separations,” Opt. Lett. 21, 546–548 (1996). [CrossRef] [PubMed]
40. G. Videen, R. G. Pinnick, D. Ngo, Q. Fu, and P. Chýlek , “Asymmetry parameter and aggregate particles” Appl. Opt. 37, 1104–1109 (1998). [CrossRef]
41. Y. -L. Xu, “Scattering Mueller matrix of an ensemble of variously shaped small particles,” J. Opt. Soc. Am. A 20, 2093–2105 (2003). [CrossRef]
42. D. Passos, J. C. Hebden, P. N. Pinto, and R. Guerra, “Tissue phantom for optical diagnostics based on a suspension of microspheres with a fractal size distribution,” J. Biomed. Opt. 10, 064036–1–11 (2005). [CrossRef]