1.

Introduction

A hospital laboratory on average examines about 300 blood samples in vitro for diagnostic purposes daily.¹ As a prominent component of blood, mature human red blood cells (RBCs) are light scatterers with homogeneous bodies enclosed by membranes and have attracted significant attention for optical diagnosis of disorders related to RBCs and blood. Flow cytometry is an established technique for automatic measurement and/or sorting of RBCs through detection of scattered light featured with good SNRs and detection simplicity and no necessity for cell staining.¹ In ectacytometry it has been demonstrated experimentally^{2 3} that the distribution of the forward light scattering depends on the size, shape, and composition of individual RBCs. To better understand existing experimental data, develop new optical diagnostic tools, and provide important information for understanding light-tissue interaction in vivo and in vitro, it is essential to quantitatively model light scattering by a single RBC.

RBCs possess intriguing biconcave shapes, which vary with the pressure and physiological and pathological conditions found in blood flow. For example, it has been reported that the biconcave shape can change toward a cup (stomatocytic) shape or a spiculated (echinocytic) shape as a result of variation in the ionic strength and pH values of the immersion medium.⁴ The size parameters of RBCs are in the range of 10 to 50 for wavelengths of general interest, where the size parameter (α) is defined as (2πR)/λ with 2R as the characteristic size of the particle and λ as the light wavelength in the surrounding medium. No analytical solutions can be found for light scattering from scatterers with these properties. Early investigations of light scattering by RBCs have used Mie theory, with the sphere as either an approximate model⁵ or an accurate model of spherical RBC obtained through chemical treatment.⁶ In ektacytometry based studies, light scattering by RBCs has been investigated with an anomalous diffraction approximation by assuming ellipsoidal shapes.⁷ Light scattering by an RBC of undeformed biconcave shape has been analyzed using a straight-ray approximation.⁸ This approximation is valid only in the small scattering angle region under the condition of small refractive index of RBC relative to the surrounding medium. Recently, a boundary-element method was developed to solve the scattering problem of axisymmetry for an undeformed RBC of biconcave shape aligned with the direction of incident light.^{9 10} The boundary-element method has the advantages of high accuracy and low computational cost for scatterers of axisymmetric shapes and homogeneous bodies. But these advantages no longer hold when it comes to general cases of light scattering from biological cells, where arbitrary shape and inhomogeneous dielectric response within the cell are common.

Cellular scattering of light at large angles toward the backscattering direction was demonstrated experimentally to contain rich information on cellular structure.¹¹ Therefore, it is highly desirable to develop reliable methods to model light scattering by individual RBCs with an arbitrary shape in the full range of scattering angles. In general, numerical methods must be used to model light scattering processes by biological cells with complex shapes and inhomogeneous bodies. Recently, Mishchenko et al.¹² comprehensively reviewed the light scattering computational methods, which include the T-matrix method,¹³ the discrete dipole approximation,^{14 15} and the finite-difference time-domain (FDTD) technique.^{16 17 18 19 20} In this paper, we use the FDTD method to compute the scattering properties of a single RBC. Although computationally intensive for three-dimensional (3-D) problems, this method employs an algorithm that is simple to adapt to a wide variety of biological cells for calculating scattered electromagnetic fields. Furthermore, the computational time of the FDTD calculations can be significantly reduced using parallel computing techniques; making it attractive as the computing power and storage size increase rapidly at a constant cost. The FDTD method has been applied to light scattering from epithelial cells to investigate the effect of intracellular components^{21 22} and to study the effects of shape change in RBCs,²³ and the results of these studies demonstrate the potential of the application of the FDTD method in cell scattering studies. As a prelude to our long-term effort in accurate modeling of light scattering from biological cells,²⁴ we investigated numerically the light scattering properties of a biconcave RBC with a serial FDTD code. In this paper, we present results on Mueller scattering matrix calculations of scattered light by an RBC of biconcave shape and their dependence on wavelength of light and the shape and orientation of the cell. The FDTD method and modeling of RBC are described in the next section, followed by results and discussions.

2.

Mathematical Modeling

2.2.

RBC Modeling

A mature human RBC is a nonnucleated cell filled with an aqueous solution of hemoglobin proteins within the membrane and is modeled in this report as a homogenous biconcave body. The hemoglobin proteins have different derivatives according to their binding molecules and we consider only the two major components of oxyhemoglobin (HbO ₂) and deoxyhemoglobin (Hb). For a homogeneous RBC, the refractive index n=n_r+in_i is treated as a constant within the cell with the imaginary part given by

Eq. (6)

$$n_i=\frac{\text{ln 10}}{4\pi }\frac{c{\lambda }_0}{n_r}\left(h_1{\gamma }_1+h_2{\gamma }_2\right),$$

where λ₀ is the light wavelength in free space, c is the hemoglobin concentration, and h₁ and h₂ (=1−h₁) and γ₁ and γ₂ are the volume ratio and molar extinction coefficients of HbO ₂ and Hb, respectively. The spectroscopy data of γ₁ and γ₂ were reported for λ₀ ranging from 200 to 1000 nm in Ref. 34. We assumed a hemoglobin oxygen saturation of h₁=97% for a normal RBC and a constant value of n_r=1.400 because the wavelength dependence of the real refractive index is not known but expected to be weak in this spectral region.⁷ In addition, we adopted c=335 g/l or 5.04×10⁻³ mole/liter for the hemoglobin concentration.³⁵ With these data, the imaginary refractive index of a normal RBC was calculated as a function of λ₀ between 400 and 1000 nm and is plotted in Fig. 1. Note here that the value of n_i calculated in this paper is about one order of magnitude smaller than that reported in Ref. 8. In calculating the scattered fields, we assumed that the RBC is suspended in blood plasma which is transparent to the light wavelength considered in this study with a constant refractive index n_p=1.350.

Figure 1

Imaginary refractive indices of Hb, HbO ₂, and a normal RBC as a function of wavelength λ₀ with h₁=97%.

Human RBCs possess a viscoelastic structure and can change shape isovolumetrically in response to different shear or physiologic conditions, resulting in significant pathological consequences.^{4 36} In modeling the shape of the RBC, we employed a simple surface function introduced in Ref. 8 to describe the biconcave disc shape of the RBC as

Eq. (7)

$$r({\theta }_A,{\varphi }_A)=a{\text{sin}}^q{\theta }_A+b,$$

where θ_A and ϕ_A are the polar and azimuthal angles in spherical coordinates relative to the axis of rotational symmetry, shown in Fig. 2(a) as the A axis, and the radius of the disk is determined by R=a+b, center thickness by 2b, and the shape by q with q=0 corresponding to a sphere. We adopted the following shape parameters: q=5, a=3.0 μm, and b=0.75 μm for a typical undeformed biconcave shape in our calculations. The RBC possessing these parameters has a volume of V=78.4 μm ³, equivalent to the volume of a sphere with a radius of R=2.66 μm. Two additional biconcave shapes with values of q=2 and 9 were also used in this report to study the effect of shape variation and the isovolumetric shape changes of an RBC was realized by adjusting a and b. The corresponding parameters for these two shapes are q=2, a=2.1 μm, b=1.12 μm, and q=9, a=3.8 μm, b=0.41 μm. The cross sections of these cases are shown in Fig. 2(b). The geometrical cross section σ_g toward the incident light was calculated numerically for each shape and orientation in the efficiency factor calculations. This was done by projecting the illuminated portion of the RBC surface (i.e., shadowed surface region excluded) onto the plane that is perpendicular to the incident light. Although the shape function given by Eq. (4) does not fit precisely to the shape of RBC determined experimentally,³⁷ it was easy to be implemented in the code and sufficiently accurate for our purpose of investigating light scattering dependence on shape, orientation and wavelength.

Figure 2

(a) Configuration of light scattering by an RBC, where the scattering plane is defined by the propagation direction of the incident light of z axis and that of the scattered light and marked by the two dashed lines; and (b) the cross sections of undeformed (q=5) and deformed RBC of biconcave disk shapes in a polar plot.

3.

Results

The scattering configuration is depicted in Fig. 2(a) with an RBC centered at the origin of the coordinate system of (x,y,z) and the incident light field propagated along the z axis. As discussed earlier, for the amplitude matrix, two incident waves with independent polarization are required. In these computations, the FDTD simulations were carried out for two incident polarization configurations with linear polarization either along the x or y axis to obtain the full amplitude matrix and subsequently the scattering matrix. The orientation of the RBC is indicated by the angle of incidence θ_i defined as the angle between the z axis and the A axis. The Mueller scattering matrix is, in general, a function of both scattering angle θ_s and azimuthal angle ϕ_s. In the case of normal incidence with θ_i=0 deg, the Mueller scattering matrix depends only on the scattering angle θ_s. We chose the same grid size along the x, y, and z directions: Δx=Δy=Δz.

The FDTD code was first validated by comparing with Mie theory for the case of a homogenous sphere suspended in blood plasma of refractive index n_p. The sphere was assumed to possess the same optical properties of a normal RBC with n_r/n_p=1.0385 from Ref. 8 and n_i=1.6804×10⁻⁵ according to Fig. 1 at the wavelength of λ₀=1 μm. The radius of the sphere R was set at 2.655 μm to obtain the same volume as an undeformed RBC. The FDTD simulation was carried out with a grid size of λ/Δx=30, where λ=λ₀/n_p is the wavelength in the blood plasma. The angular distribution of S₁₁, or the single-scattering phase function, of the Mueller scattering matrix calculated by the FDTD method is compared with that of Mie theory in Fig. 3 in a semilog scale. Errors in the FDTD calculation in comparing with Mie theory for the four nonzero independent elements of the Mueller scattering matrix are presented in Fig. 4, where the error of S₁₁ is in the form of relative error and that of other elements is in the form of absolute error. The relative error of S₁₁ is defined as the difference of calculated S₁₁ between the two methods scaled by that of Mie theory while the absolute errors of the other elements are the differences between the results of the two methods. In general, deviations of the FDTD results from Mie theory become evident for three scattering angles large than 90 deg that are associated with local minima in the phase function. The largest error is observed at θ_s=151 deg; where the scattered energy is at the lowest level. The relative errors are not presented for elements other than S₁₁ for these elements have values of 0 or close to 0 at some scattering angles and thus relative errors cannot be defined.

Figure 3

Mueller scattering matrix element S₁₁ calculated by FDTD and Mie theory for a homogenous sphere with radius R=2.655 μm, n_r/n_p=1.0385, and n_i=1.6804×10⁻⁵ at wavelength λ₀=1 μm. Grid size factor for the FDTD calculation is given by Δr=λ/30.

Figure 4

Errors in the four independent Mueller scattering matrix elements calculated by the FDTD method compared with the Mie theory for the same calculation in Fig. 3.

The computational resource requirements (memory and CPU time) for the 3-D FDTD problem increase with grid size according to a power law. It is observed in the current calculation that the CPU time increases approximately with the third power of grid size while the memory requirements increase approximately with a power of 2.4. The choice of a suitable grid size is important to improve the code performance. To evaluate the effects of grid size, we calculated the Mueller scattering matrix of an undeformed RBC with two different grid sizes of λ/Δx=30 and λ/Δx=20. The results demonstrate that the effect of grid size change from λ/Δx=30 to λ/Δx=20 on the calculated angular distributions of light scattered by an undeformed RBC at θ_i=0 deg is negligible except near θ_s=113 deg and θ_s=176 deg, where the maximum relative errors near these two angles (±4 deg) were found to be 130 and 150%, respectively. However, with a grid size of λ/Δx=20, the CPU time and memory use can be reduced by a factor of 3 in comparison with λ/Δx=30. In addition, the interest of current study is on the dependence of angular distribution of light scattering on the RBC shape, orientation and light wavelength, simulations with λ/Δx=20 provides sufficiently accurate results for this aspect. For these reasons, we used the grid size of λ/Δx=20 for the rest of the calculations presented in this paper.

We first applied the FDTD method to investigate the effects on the angular distributions of scattered light and optical parameters caused by isovolumetric shape changes of a RBC by varying the shape parameter q in Eq. (4) from the undefromed case of q=5. The Mueller scattering matrix elements have been calculated using the FDTD method for each of the three biconcave shapes in Fig. 2(b) for a normally incident light of θ_i=0 with the same real and imaginary indices of refraction used for the sphere case at λ₀=1 μm. It is observed that the shape change correlates with scattering pattern change, and marked changes are evident in the phase functions S₁₁ and S₄₃. As the shape deviates from the spherical shape by isovolumetrically “flattening” along the z axis, the phase function shows drastically enhanced backscattering near θ_s=180 deg by factors of more than 10, while side scattering near θ_s=90 deg demonstrates large pattern shifting. The results for q=2 and 9 are plotted in Fig. 5 for the four independent elements of the Mueller scattering matrix to illustrate the pattern. The corresponding optical parameters are compared in Table 1 together with those of the sphere case of q=0, and it shows large drops in scattering and absorption power as the RBC is flattened for increasing q, but the correlation between g and q is relatively weak.

Figure 5

Mueller scattering matrix elements as a function of scattering angle θ_s at shape parameters of q=2 and 9 for a biconcave RBC with θ_i=0 deg and λ₀=1 μm.

Table 1

We next studied the effect of RBC’s orientation on the Mueller scattering matrix elements by considering nonzero values of incident angle θ_i. In these cases, the system loses its symmetry with respect to the z axis, and the Mueller scattering matrix elements depend on both θ_s and ϕ_s. In addition, the number of independent elements of the Mueller scattering matrix can increase to 6. For this part of the calculations, we considered only the undeformed shape of RBC with q=5 at λ₀=1 μm. The phase functions S₁₁(θ_s,ϕ_s) for different incidence angles θ_i of 5, 45, and 90 deg are depicted in Fig. 6. Here we observe a clear dependence of the scattering pattern on ϕ_s as the A axis tilted away from the z axis [see Fig. 2(a)] and a return to increased symmetry at θ_i=90 deg. Particularly, a relatively large portion of the energy is scattered near ϕ_s=90 deg for the case of θ_i=45 deg. The azimuthally averaged Mueller scattering matrix was obtained and four of the independent elements for different incident angles are plotted in Fig. 7 as a function of the scattering angle θ_s. We can see that the averaging over the azimuthal angle ϕ_s smoothes the angular oscillations in the elements for nonzero θ_i comparing with the case of θ_i=0 deg. The differences between the element S₂₂ and S₁₁ and S₄₄ and S₃₃ were found to be very small and thus are not presented in Fig. 7. The dependency of the optical parameters on the incident angle is plotted in Fig. 8, and it shows a clear pattern of increasing in scattering and absorption powers as the angle of incident increases, and the decrease in g indicates that more energy is scattered into scattering angles other than the forward direction. This effect is demonstrated in the phase function by the broadened peak in the small scattering angle region for larger angles of incidence.

Figure 6

Contour plots of the phase function S₁₁ as a function of θ_s and ϕ_s at different values of incident angle θ_i for a typical RBC with q=5 and λ₀=1 μm.

Figure 7

Mueller scattering matrix elements as functions of the scattering angle θ_s for different angles of incidence θ_i together with the corresponding elements averaged over the incident angle for a typical RBC with q=5 and λ₀=1 μm. The elements were averaged over the azimuthal angle ϕ_s for cases of nonzero θ_i. The Henyey-Greenstein phase function calculated using the g averaged over the incident angle is presented with S₁₁.

Figure 8

Scattering efficiency Q_s, absorption efficiency Q_a, and anisotropy factor g as functions of the incident angle θ_i for the case defined in Fig. 7. The solid lines are meant only as guides for the eye.

For biconcave RBCs suspended in blood plasma, the orientations of the cells may be arbitrary. Experimental measurements of single scattering light distribution from a cell suspension correspond in general to calculated ones averaged over ϕ_s and θ_i. For this reason, we averaged the azimuthally averaged elements of the Mueller scattering matrix over eight incident angles between 0 and 90 deg, which resulted in smoothed curves that are also plotted in Fig. 7. The anisotropy factor corresponding to the averaged phase function is g¯=0.9888. The averaged phase function was calculated by weighting the azimuthally averaged phase function for each incident angle with its corresponding scattering cross section before normalization according to Eq. (5). The averaged efficiencies were also computed by obtaining the ratio of averaged cross sections to the averaged geometric cross section, and it yielded the average values of Q¯_s=0.8091 and Q¯_a=8.36×10⁻⁴.

In the study of multiple light scattering in turbid tissues, a single-parameter function of Henyey-Greenstein (HG) is widely used as the phase function for single scattering:

We also examined the scattering properties of an undeformed RBC with q=5 at different values of wavelength λ₀. It is well known that scattering and absorption cross sections are closely related to the complex index of refraction and that the frequency of oscillation in the phase function depends on the size parameter α. To study the effect of the wavelength on light scattering by an RBC, we chose two additional wavelengths corresponding to a maximum and a minimum in the imaginary index n_i between 400 and 1000 nm (see Fig. 1): λ₀=0.575 μm with n_i=5.4×10⁻⁴ and λ₀=0.700 μm with n_i=4.3×10⁻⁶. The real indices are taken to be n_r/n_p=1.0386 and 1.0385, respectively, according to Ref. 8. The Mueller scattering matrix elements for a typically shaped RBC of q=5 with a normal incidence of θ_i=0 deg at these two wavelengths are calculated and the results for λ₀=0.575 μm is compared in Fig. 9 with the case of λ₀=1 μm to show the general trend while the optical parameters are listed in Table 2 together with those of λ₀=1 μm. It is well known from Mie theory that in the case spheres, an increasing size parameter α leads to more oscillation in the angular distribution of phase function. This trend is confirmed by the results in Fig. 9 and those of the λ₀=0.7 μm case. These results also demonstrate a previously known pattern that at shorter wavelengths (i.e., large size parameter), the RBC scatters more light near the forward direction and thus decreases the scattering at large angles, causing a larger g. The optical parameters listed in Table 2 show that the scattering efficiency of an RBC increases with decreasing wavelength, and the absorption efficiency of an RBC relates positively with the imaginary part of the index (n_i) with the smallest value at λ=0.700 μm and the largest at λ=0.575 μm.

Figure 9

Comparisons of the angular distribution of the Mueller scattering matrix elements at wavelengths λ₀=0.575 μm and 1 μm for the case of θ_i=0 deg and q=5.

Table 2

4.

Discussion

In this paper, we validated a serial FDTD code for accurate modeling of light scattering by a biconcave RBC and obtained the Mueller scattering matrix elements as a function of scattering angles. Our results reveal a direction correlation between angular distribution of Mueller scattering matrix elements and the shape and orientation of a biconcave RBC at different light wavelengths. It was shown that the scattering and absorption efficiencies decrease as the RBC is flattened, and those parameters increase as the incident light tilts from the axis of symmetry of the cell. Results obtained at different wavelengths indicate that the scattering efficiency increases rapidly with decreasing wavelength, while the absorption efficiency depends sensitively on the imaginary part of the refractive index of the cell.

As presented in Figs. 3 and 4 for the case of a sphere, the errors in the Mueller scattering matrix calculations in our FDTD method are significant only at a few isolated scattering angles where the scattered fields are very weak, and the deviation of FDTD results from the Mie theory is often dominated by errors near the backscattering direction of θ_s=180 deg. Although similar errors are also observed in other calculations,^{17 21} it appears more pronounced in the large scattering angle region for scatterers like biological cells in a fluid, particularly near the backscattering direction. This may be attributed to the fact that the relative index of refraction in the case of light scattering from biological cells in a fluid is very close to 1. In this case, the phase function is featured with a strong peak in the forward direction, which drops rapidly with increasing scattering angle and a dynamic range of more than seven orders of magnitude. This feature of the phase function subjects the low-scattering-energy region, namely, the large scattering region, to higher relative errors. The sources of errors in this calculation may come from the use of large grid size Δx, shorted time marching duration ΔT, and deficiency in the absorbing boundary condition. Large Δx causes imperfect fitting of the spherical surface with the finite sized grid used in the calculations, while small ΔT leaves some measurable scattered energy uncounted.¹⁷ Reducing either grid size or increasing time-marching duration can increase the accuracy of FDTD simulations at the expense of higher computational cost. The PML absorption boundary condition used in this calculation is very efficient for most calculations, but for the reason discussed earlier, it seems not sufficient for the large dynamic range of the phase function in this case. A recently introduced uniaxial perfectly matched layer (UPML) absorption boundary condition^{16 38} reportedly may provide improved efficiency for certain cases. We are currently testing this possibility by implementing UPML into our program.

Since a mature human RBC contains no internal structure, the angular distribution of S₃₄ does not display systematic patterns of changes for different shapes, orientations, and wavelengths, as observed in some cases where internal structure was presented.^{29 33} On the other hand, it is interesting to notice that the phase function S₁₁ in Fig. 5 exhibits a very strong enhancement near θ_s=180 deg when the shape of the biconcave RBC deviates from the typical undeformed shape of q=5. As discussed earlier, the use of a relatively large grid size of λ/Δx=20 may contribute up to 150% simulation errors in this region, or a factor of 1.5, the changes in the backscattering probabilities, however, are much larger than that when q changes either to 2 or 9. If the S₁₁ is integrated from θ_s=170 to 180 deg to obtain (S₁₁)_b, we found that (S₁₁)_b of an RBC increases from 0.0712 for the case of q=5 to 1.14 for q=2 and 4.35 for q=9. Furthermore, one can observe from Fig. 7 that the S₁₁ for the case of q=5 changes very little near θ_s=180 deg for different angle of incidence θ_i. These results in combination suggest that enhanced backscattering from RBCs of biconcave shapes may be used as a sensitive tool to measure changes in RBC shape.