Open Access
Add to BibSonomyAbstract
The discrete sources method (DSM) and the discrete dipole approximation (DDA) were compared for simulation of light scattering by a red blood cell (RBC) model. We considered RBCs with diameters up to 8 μm (size parameter up to 38), relative refractive indices 1.03 and 1.06, and two different orientations. The agreement in the angle-resolved S 11 element of the Mueller matrix obtained by these methods is generally good, but it deteriorates with increasing scattering angle, diameter and refractive index of a RBC. Based on the DDA simulations with very fine discretization (up to 93 dipoles per wavelength) for a single RBC, we attributed most of the disagreement to the DSM, which results contain high-frequency ripples. For a single orientation of a RBC the DDA is comparable to or faster than the DSM. However, the relation is reversed when a set of particle orientations need to be simulated at once. Moreover, the DSM requires about an order of magnitude less computer memory. At present, application of the DSM for massive calculation of light scattering patterns of RBCs is hampered by its limitations in size parameter of a RBC due to the high number of harmonics used for calculations.
©2010 Optical Society of America
1. Introduction
There is a growing interest in human blood cells, due to an important role they play in human health. Biophysical properties of blood cells are sensitive markers for different diseases [1]. Some of the efficient modern methods of blood cells analysis are based on light-scattering measurements. Spatial distribution and polarizing properties of light scattered by a cell depend on its morphology: shape, size, refractive index, and internal structure [2]. Blood cell analysis is based on a solution of the inverse light-scattering problem (ILSP), so the accuracy of the latter determines efficiency of the former. Mathematical modeling of light scattering from a cell, i.e. solution of the direct problem, is an essential stage in development of a solution of an ILSP. Since ILSP for a single particle does not have a closed form solution, one has to numerically invert a mapping from parameter space to space of light scattering signals, based on the approximate description of this mapping. For the latter the direct problem has to be solved multiple times (e.g. thousands), each requiring a significant computer resources, because blood cells are much larger than the wavelength of the visible light and have complex morphology. Therefore, comparison of different codes for simulation of light scattering by blood cells in terms of accuracy, computational time, and memory requirements is a relevant problem.
Several research groups worldwide have simulated light scattering by a red blood cell (RBC) [3–7]. The RBC attracts special interest because it is a homogeneous particle with a shape varying from a nearly spherical through a biconcave disk to a toroidal one. Therefore, an important problem for simulation is to find an appropriate cell model reproducing the real cell shape. It is further complicated by variation of refractive index, depending on hemoglobin concentration, and relatively large size parameter (typically in the range of 30-60, when considering visible light). However, the RBC shape has a symmetry axis, which expands a set of applicable light scattering methods. During the last years the following methods were used to simulate light scattering by a single RBC: boundary element method [3], finite difference time domain method [4,6], discrete dipole approximation (DDA) [5,6], and discrete sources method (DSM) [8–10]. Earlier approximate shape models for RBC were used, e.g. oblate spheroid [11], which allowed simpler light scattering methods to be applied. Several of these methods were compared [7] and it was concluded that DSM is the most efficient method, combining the fastest computational speed with the ability to calculate light scattering for any directions of incident light at once. However, the comparison was based on a RBC with size parameter 28. Moreover, the light was incident along the RBC symmetry axis, which led to a considerable performance gain for the DSM.
We do not consider here approximate light scattering methods, e.g. Wentzel-Kramers-Brillouin and Rytov approximations [6,12], since they produce satisfactory results only for particular light scattering quantities (e.g. near-forward scattering). The methods described above are rigorous in the sense that potentially they can be made as accurate as required given infinite computer power. However, they do have practical limitations, due to either enormous requirements for memory and computational time or numerical instability for a fixed arithmetic precision.
In this manuscript we perform a direct comparison of two most promising methods (DDA and DSM) in simulation of light scattering from a single RBC modeled by a biconcave disk with size parameter up to 38. We compare the angular dependency of the element S 11 of the Mueller matrix calculated with both methods varying cellular characteristics within typical biological ranges. A special attention is paid to limitations of both methods in terms of size parameter of a RBC.
2. Simulation methods
2.1 The discrete sources method
The DSM theory and its numerical scheme has presented in previous papers [9,13], here we are just shortly discussing some aspects. The DSM uses an idea of approximate solution which is constructed by representing the electromagnetic fields as a finite linear combination of the electric and magnetic fields of multipoles distributed inside the scatterer. In case of RBC discrete sources (DS) are deposited in a complex plane adjoined to the symmetry axis of the particle. This procedure is presented in [14] in details.
Due to the special construction of the approximate solution, which takes into account the rotational symmetry of the scatterer and deposition of DS in the complex plane adjoined to the symmetry axis, the matching of the transmission conditions at the particle surface is reduced to a set of one-dimensional problems enforced at the particle generatrix. It has been found that more stable results can be obtained by using the generalized point-matching technique and a pseudo-solution of an over-determined system of linear equations [14].
DSM provides the opportunity to control the accuracy of the computational result within two steps: (a) сontrol of the internal convergence of the results by increasing the number of matching points and DS and (b) checking the residual in least square norm of the boundary conditions at the particle surface. In a simulation run the number of matching points where the DS amplitudes are defined can be increased until the maximum accuracy of the results is achieved. The DS number usually is 2-4 times less then the number of matching points. The detailed numerical scheme for choosing matching points is described in [14].
2.2 The discrete dipole approximation
The DDA is based on volume discretization of the particle and thus is applicable to scatterers of arbitrary shape and internal structure. A couple of extensive reviews of the method are available [15,16]. The main limitation of the DDA is its computational complexity, which grows proportional to the number of volume discretization elements (dipoles). As a numerical implementation of the DDA we have used ADDA v.0.79 [17], which is capable of running on a cluster of computers parallelizing a single DDA computation [18]. We used the default discretization in ADDA, corresponding to 10-11 dipoles per wavelength for RBCs considered in this manuscript. Being a general method, the DDA do not employ particle symmetries except one particular case. When particle is invariant under rotation by 90° over the propagation direction of the incident light, the complete Mueller matrix can be obtained from simulation for only one incident polarization [18]. Then, simulation time is twice less than for a general case.
3. Optical model of a RBC
A mature RBC has a biconcave discoid shell composed by hemoglobin (32%), water (65%), and membrane components (3%), all in mass percentage, and does not contain any nucleus [19]. The shape of the RBC can be described by several formulae introduced by Skalak et al. [20], Fung et al. [21] and Yurkin [22, Chapter 4.3]. In the current study the DSM and the DDA were applied to simulate light scattering from an individual RBC with the widely used cell shape model introduced by Fung et al. [21]:
where z and ρ are cylindrical coordinates and D = 2R is the diameter of a RBC (the only free parameter in this model). A particular case of D = 6 μm is shown in Fig. 1 . The morphology data for RBCs from different sources has been reviewed by Yurkin [22, Chapter 4.1]. The diameter of a RBC typically varies from 6 to 9 μm, and the real part of the refractive index at λ = 0.66 μm falls between 1.39 and 1.41. Imaginary part of the refractive index is about 10−4 [23], therefore we neglect it in the current study. We also do not consider RBC outer membrane, which has a thickness of 7 nm [24].We compare the simulation methods using the following characteristics of a RBC. Diameter is varied from 6 to 8 mm, the upper limit being determined by the current capabilities of the DSM (discussed in Section 4). We consider RBCs suspended in buffered saline (refractive index 1.337). Then two considered values of the relative refractive index m are 1.03 and 1.06, and wavelength in the medium is fixed at λ = 0.4936 μm (0.66 mm in vacuum). Size parameter of the simulated RBCs is from 28 to 38. Incident radiation propagates along the z-axis, and we calculate dependence of S 11 on the scattering angle θ in the xz-plane. Orientation of a RBC is with symmetry axis along either z-axis or x-axis (the Euler angle β equals 0° or 90° respectively).
Since the DDA and the DSM are volume- and surface-based methods respectively, they work with different descriptions of the particle model. To make our results independent from the shape description we have started from the same RBC profile described by 100 nodes, assuming linear interpolation between the nodes. The DSM uses this profile as is, while the DDA constructs volume discretization of a RBC from this profile. The latter is automatically done by ADDA using the “-shape axisymmetric” command line option [25].
4. Results and discussion
Figures Fig. 2 -Fig. 5 present S 11(θ) calculated using the DSM and the DDA varying diameter, refractive index and orientation of a RBC. The agreement between the two methods is generally good however it strongly depends on problem parameters. General tendency is that differences increases with θ, D, β, and m. Refractive index has the least effect on the difference, when varied inside biological range for a RBC. For β=0° the agreement is good along the whole θ range when D ≤ 7 μm, but there is significant (order of magnitude) disagreement for 90° ≤ θ ≤ 120° when D = 7.5 μm. For β=90° and D ≤ 7.5 μm the agreement is good only up to θ=60°-70°. The largest tested diameter (8 μm) shows principal disagreement for all test cases, except the near-forward scattering (up to 15).
Fig. 2 S 11(θ) in logarithmic scale for RBCs with D = 6 mm, m = 1.03 (a,c) and 1.06 (b,d), β = 0° (a,b) and 90° (c,d) simulated with the DSM and the DDA.
Since we compare two potentially inaccurate methods, it is desirable to have some reference solution to judge the accuracy of both methods independently. However, to the best of our knowledge no principally more accurate methods exist for a RBC. Moreover, the accuracy of the DSM cannot be easily improved by using extra computational resources. For the oblate scatterers of large diameter DSM scheme requires calculation of many Fourier harmonics. For example, for RBC with diameter of 6 μm and exciting wavelength λ = 0.4936 μm 45 harmonics are needed. This causes numerical instability for large orientation angles β. For example, the surface residual for all RBC diameters presented in paper does not exceed 0,2% for β=0°. At the same time the surface residual for β = 90° varies from 3% for D = 6 µm to 8% for D = 8 µm. Although the surface residual can be used as an internal quality test, it is only an approximate measure of the accuracy of the final scattering quantities.
Therefore, we invested additional computational resources in improving the accuracy of the DDA, relying on the proven convergence of the DDA results with refining discretization [26]. We performed DDA simulations for a single RBC (D = 7.5 μm, m = 1.03) and two orientations (β = 0° and 90°) varying dipole size d from λ/8 to λ/93. Number of dipoles per grid varied from 128 to 1408 respectively, and total number of dipoles was up to 6×108. These huge simulations were carried out on the Dutch compute cluster LISA [27] using up to 560 processor cores and up to 770 GB of memory. The results for a particular scattering angle (θ= 120°) are shown in Fig. 6 . Although this is one of the worst convergence among all θ (data not shown), one can clearly see that the DDA results indeed converge with decreasing dipole size. However, this convergence is oscillating, complicating the choice of a particular reference value or interval. For instance, using extrapolation technique as described by Yurkin et al. [28], based on 9 best discretizations, leads to confidence interval [0.54,3.14]×10−2 for S 11(120°) for the case of β = 0°, which is too wide for practical purposes. The extrapolation technique was originally tested on scatterers with wavelength-sized scatterers [28] and it requires a separate study to be effectively applied to much larger scatterers.
Fig. 6 Convergence of the DDA results for S 11(120°) with refining discretization for a RBC with D = 7.5 mm, m = 1.03, β = 0° (a) and β = 90° (b). Confidence intervals determined from these data (see text) are shown by hatched bands.
In this manuscript we use a simpler approach – confidence interval is determined by maximum and minimum of DDA results for dipole sizes less than λ/16. Such confidence intervals for S 11(120°) are [1.59,2.43]×10−2 and [2.06,2.32]×10−1 for β = 0° and 90°respectively (also shown in Fig. 6). One can see that they seem reliable, i.e. further decrease of dipole size should not move the DDA results out of these intervals. Similar conclusions can be reached for other θ (data not shown), however, this reliability is empirical rather than rigorously proven. The confidence bounds obtained by this method for all scattering angles are shown in Fig. 7 together with DSM and DDA (with default discretization) results for the same RBCs. Confidence bounds are very narrow and coincident the DDA results for default discretization over the major part of the angular range. Assuming the reliability of the confidence bounds, we conclude that the DSM is significantly less accurate than the DDA for this particular RBC, at least in the range of θ with the largest difference between the two methods. Accuracy of the DDA is expected to have small dependence on particle size for fixed dipole size [18]. Therefore, we conclude that most of the difference between the results of the two methods can be attributed to the inaccuracy of the DSM, in agreement with DSM internal error estimates.
Fig. 7 Same as Fig. 4 (a,c), but with addition of confidence bounds obtained by DDA simulations with dipole sizes from λ/17 to λ/93 (see text).
Next, we turn to computational considerations. All production runs, results of which are shown in Figs. 2–5, were run on the same computer – AMD Athlon 64 × 2 Dual Core 3800 + , 2.01 GHz with 2 GB RAM. Both codes were run in sequential mode (on a single core). Performance results are presented in Table 1 . To make a fare comparison we covered a complete angle in one scattering plane with step of 0.5° (totally 720 scattering angles) with both methods. Additionally to two RBC orientations we consider simultaneous calculation for 10 values of β (from 0° to 90° with step 10°), which is a typical task when constructing a database of light scattering patterns [5,22]. The DDA time for the latter case were estimated from measured times t(0°) and t(90°) for two orientations as 9t(0°) + 5t(90°), assuming linear dependence of the number of iterations in the DDA on β. The latter implies that simulation time for β ≠ 0° linearly changes from 2t(0°) to t(90°) when β increases from 0° to 90°. Memory usage for the DSM method is approximately 60 MB for all studied cases (not shown in Table 1). Memory requirements of the DDA increase with D and are about an order of magnitude larger than that of the DSM. However, it is still small enough to fit into a standard desktop computer.
Table 1. Time and memory requirements for DSM and DDA methods.a
The dependence of computational time on m is markedly different for the compared methods. The DSM speed is comparable or even slightly faster for m = 1.06 than for 1.03, while the DDA is about 50% slower for the larger m. This agrees with previous observations that the DDA performance is especially good for index-matching particles [18,29]. Apart from that the DDA and the DSM show comparable speed for β = 0°, which is a special symmetric case for both methods (see Section 2). For a single non-symmetric orientation (β = 90°) the DDA is 1.4-3.4 times faster. However, the DSM is 1.5-3 times faster for the set of 10 orientations.
The speed of the DDA, and ADDA in particular, can be improved by running it in a parallel mode, employing all available processor cores. About 50% increase in speed is expected using two cores [25], but it was not tested in this manuscript. Moreover, simultaneous DDA runs for different particle orientations can be accelerated using the ideas proposed by Okada [30], but they are not yet implemented in any publicly available DDA code. The performance of the DSM for larger oblate scatterers can be improved by use of iterative solvers for matrix pseudo-inversion and optimization of DS distribution.
5. Conclusion
The DSM and the DDA were compared for simulation of light scattering by a RBC model with diameter up to 8 μm (size parameter up to 38). The most surprising result is that overall the DDA successfully competes with if not outperforms the DSM, in spite of the axisymmetric shape of a RBC. That is probably due to the lower relative refractive index of (from 1.03 to 1.06) which results fast convergence of the iterative scheme. For such “weak” scatterers even Born approximation provides reasonable results for near-forward scattering direction.
The agreement in the angle-resolved S 11 element of the Mueller matrix between the two methods is generally good, however it deteriorates with increasing scattering angle, diameter and refractive index of a RBC, and when switching from symmetric to non-symmetric orientation of a RBC with respect to the incident radiation. Separate convergence study of the DDA for a single RBC involving up to 6 × 108 dipoles showed that most of the disagreement between the methods can be attributed to the DSM. The relatively worse accuracy of the DSM is also evidenced by its internal error estimates and high-frequency ripples present in its results, especially for larger RBC diameter.
For a single orientation of a RBC the DDA is comparable to or faster than the DSM. The important advantages of the DSM are that it requires about an order of magnitude less computer memory and that it is faster when a set of particle orientations are considered at once. The latter property is favorable for construction of a database of light scattering patterns for different sizes, refractive indices, and orientations of a RBC, e.g. to solve an inverse light scattering problem. However, application of the DSM to this task is currently hampered by its limitations in size parameter of a RBC (up to about 40) due to numerical instability. In contrast, the DDA, in particular ADDA code, is easily applicable to larger RBC if sufficient memory is available (either shared or distributed among several computers).
Acknowledgements
This research is supported by grants of the Deutsche Forschungsgemeinschaft (DFG)ER553/1-1, Russian Foundation for Basic Research, No 07-04-00356-a, and No 08-02-91954-NNIO_а, integration grants of the Siberian Branch of the Russian Academy of Science, No 2009-37, and No 2009-7, grant from the program of Presidium of the Russian Academy of Science, No 2009-27b, and by program of the Russian Government “Research and educational personnel of innovative Russia” (contract P2497).
References and links
1. J. P. Greer, J. Foerster, and J. N. Lukens, eds., Wintrobe's Clinical Hematology, (Lippincott Williams & Wilkins Publishers, Baltimore, USA, 2003).
2. V. P. Maltsev, and K. A. Semyanov, Characterisation of Bio-Particles from Light Scattering, Inverse and Ill-Posed Problems Series (VSP, Utrecht, 2004).
3. S. V. Tsinopoulos and D. Polyzos, “Scattering of He-Ne laser light by an average-sized red blood cell,” Appl. Opt. 38(25), 5499–5510 (1999). [CrossRef]
4. J. Q. Lu, P. Yang, and X. H. Hu, “Simulations of light scattering from a biconcave red blood cell using the finite-difference time-domain method,” J. Biomed. Opt. 10(2), 024022 (2005). [CrossRef] [PubMed]
5. M. A. Yurkin, K. A. Semyanov, P. A. Tarasov, A. V. Chernyshev, A. G. Hoekstra, and V. P. Maltsev, “Experimental and theoretical study of light scattering by individual mature red blood cells by use of scanning flow cytometry and a discrete dipole approximation,” Appl. Opt. 44(25), 5249–5256 (2005). [CrossRef] [PubMed]
6. A. Karlsson, J. P. He, J. Swartling, and S. Andersson-Engels, “Numerical simulations of light scattering by red blood cells,” IEEE Trans. Biomed. Eng. 52(1), 13–18 (2005). [CrossRef] [PubMed]
7. T. Wriedt, J. Hellmers, E. Eremina, and R. Schuh, “Light scattering by single erythrocyte: comparison of different methods,” J. Quant. Spectrosc. Radiat. Transf. 100(1-3), 444–456 (2006). [CrossRef]
8. E. Eremina, Y. Eremin, and T. Wriedt, “Analysis of light scattering by erythrocyte based on discrete sources method,” Opt. Commun. 244(1-6), 15–23 (2005). [CrossRef]
9. E. Eremina, J. Hellmers, Y. Eremin, and T. Wriedt, “Different shape models for erythrocyte: Light scattering analysis based on the discrete sources method,” J. Quant. Spectrosc. Radiat. Transf. 102(1), 3–10 (2006). [CrossRef]
10. E. Eremina, “Light scattering by an erythrocyte based on Discrete Sources Method: shape and refractive index influence,” J. Quant. Spectrosc. Radiat. Transf. 110(14-16), 1526–1534 (2009). [CrossRef]
11. A. M. K. Nilsson, P. Alsholm, A. Karlsson, and S. Andersson-Engels, “T-Matrix Computations of Light Scattering by Red Blood Cells,” Appl. Opt. 37(13), 2735–2748 (1998). [CrossRef]
12. A. N. Shvalov, J. T. Soini, A. V. Chernyshev, P. A. Tarasov, E. Soini, and V. P. Maltsev, “Light-scattering properties of individual erythrocytes,” Appl. Opt. 38(1), 230–235 (1999). [CrossRef]
13. Yu. Eremin, “The method of discrete sources in electromagnetic scattering by axially symmetric structures,” J. Commun. Technol. Electron. 45(Suppl.2), 269–280 (2000).
14. Y. Eremin, N. Orlov, and A. Sveshnikov, “Models of electromagnetic scattering problems based on discrete sources method” in: Generalizes Multipole Techniques for Electromagnetic and Light Scattering, T. Wriedt, ed. (Elsevier, Amsterdam, 1999), Chapter 4.
15. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11(4), 1491–1499 (1994). [CrossRef]
16. M. A. Yurkin and A. G. Hoekstra, “The discrete dipole approximation: An overview and recent developments,” JQSRT 106, 558–589 (2007).
17. “ADDA - light scattering simulator using the discrete dipole approximation”, http://code.google.com/p/a-dda/ (2009).
18. M. A. Yurkin, V. P. Maltsev, and A. G. Hoekstra, “The discrete dipole approximation for simulation of light scattering by particles much larger than the wavelength,” J. Quant. Spectrosc. Radiat. Transf. 106(1-3), 546–557 (2007). [CrossRef]
19. P. Mazeron, S. Muller, and H. El Azouzi, “Deformation of erythrocytes under shear: a small-angle light scattering study,” Biorheology 34(2), 99–110 (1997). [CrossRef] [PubMed]
20. R. Skalak, A. Tozeren, R. P. Zarda, and S. Chien, “Strain energy function of red blood cell membranes,” Biophys. J. 13(3), 245–264 (1973). [CrossRef] [PubMed]
21. Y. C. Fung, W. C. Tsang, and P. Patitucci, “High-resolution data on the geometry of red blood cells,” Biorheology 18(3-6), 369–385 (1981). [PubMed]
22. M. A. Yurkin, “Discrete dipole simulations of light scattering by blood cells” PhD thesis, (University of Amsterdam, 2007).
23. K. A. Semyanov, P. A. Tarasov, J. T. Soini, A. K. Petrov, and V. P. Maltsev, “Calibration-free method to determine the size and hemoglobin concentration of individual red blood cells from light scattering,” Appl. Opt. 39(31), 5884–5889 (2000). [CrossRef]
24. D. H. Tycko, M. H. Metz, E. A. Epstein, and A. Grinbaum, “Flow-cytometric light scattering measurement of red blood cell volume and hemoglobin concentration,” Appl. Opt. 24(9), 1355–1365 (1985). [CrossRef] [PubMed]
25. M. A. Yurkin, and A. G. Hoekstra, “User manual for the discrete dipole approximation code ADDA v.0.79,” http://a-dda.googlecode.com/svn/tags/rel_0_79/doc/manual.pdf (2009).
26. M. A. Yurkin, V. P. Maltsev, and A. G. Hoekstra, “Convergence of the discrete dipole approximation. I. Theoretical analysis,” J. Opt. Soc. Am. A 23(10), 2578–2591 (2006). [CrossRef]
27. “Description of the national compute cluster Lisa,” https://subtrac.sara.nl/userdoc/wiki/lisa/description (2009).
28. M. A. Yurkin, V. P. Maltsev, and A. G. Hoekstra, “Convergence of the discrete dipole approximation. II. An extrapolation technique to increase the accuracy,” J. Opt. Soc. Am. A 23(10), 2592–2601 (2006). [CrossRef]
29. M. A. Yurkin, A. G. Hoekstra, R. S. Brock, and J. Q. Lu, “Systematic comparison of the discrete dipole approximation and the finite difference time domain method for large dielectric scatterers,” Opt. Express 15(26), 17902–17911 (2007). [CrossRef] [PubMed]
30. Y. Okada, I. Mann, I. Sano, and S. Mukai, “Acceleration of the iterative solver in the discrete dipole approximation: Application to the orientation variation of irregularly shaped particles,” J. Quant. Spectrosc. Radiat. Transf. 109(8), 1461–1473 (2008). [CrossRef]
