Abstract

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 S11 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.

]. 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

]. 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

] 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.

], 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 S11 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

], 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

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

]. 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

]. 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

, 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

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 S11 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

present S11(θ) 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).

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

. 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 S11(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

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 S11(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

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).

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.

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

Table 1. Time and memory requirements for DSM and DDA methods.a

Diameter D, mm

Orientation angle β

m = 1.03

m = 1.06

DSM

DDA

DSM

DDA

Time, s

Time, s

RAM, MB

Time, s

Time, s

RAM, MB

6

0°

70

71

391

63

121

432

90°

493

144

439

290

10 angles

510

1359b

456

2539b

7

0°

92

109

639

100

205

708

90°

728

247

829

512

10 angles

749

2216b

885

4405b

7.5

0°

125

154

778

110

255

816

90°

1012

327

959

631

10 angles

1025

3021b

986

5450b

8

0°

162

184

924

147

303

997

90°

1462

443

1301

908

10 angles

1465

3871b

1344

7267b

a Memory requirements of the DDA do not depend on orientation, so only one value per RBC diameter is shown.

b This value is estimated using the results for β = 0° and 90° (see text).

. 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

]. 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.

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

]. 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 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 S11 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).

Table 1. Time and memory requirements for DSM and DDA methods.a

Diameter D, mm

Orientation angle β

m = 1.03

m = 1.06

DSM

DDA

DSM

DDA

Time, s

Time, s

RAM, MB

Time, s

Time, s

RAM, MB

6

0°

70

71

391

63

121

432

90°

493

144

439

290

10 angles

510

1359b

456

2539b

7

0°

92

109

639

100

205

708

90°

728

247

829

512

10 angles

749

2216b

885

4405b

7.5

0°

125

154

778

110

255

816

90°

1012

327

959

631

10 angles

1025

3021b

986

5450b

8

0°

162

184

924

147

303

997

90°

1462

443

1301

908

10 angles

1465

3871b

1344

7267b

a Memory requirements of the DDA do not depend on orientation, so only one value per RBC diameter is shown.

b This value is estimated using the results for β = 0° and 90° (see text).

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).

Cited By

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.