## Abstract

The goal of this work is to develop a computational framework to rapidly simulate the light scattering response of multiple red blood cells. Because the wavelength of visible light (3.8×10^{−7} m≤*λ*≤7.2×10^{−7} m) is approximately an order of magnitude smaller than the diameter of a typical red blood cell scatterer (*d*≈8×10^{−6} m), geometric ray-tracing theory is applicable, and can be used to quickly ascertain the amount of optical energy, characterized by the Poynting vector, that is reflected and absorbed by multiple red blood cells. The overall objective is to provide a straightforward approach that can be easily implemented by researchers in the field, using standard desktop computers. Three-dimensional examples are given to illustrate the approach and the results compare quite closely to experiments on blood samples conducted at the Children's Hospital Oakland Research Institute (CHORI).

## 1. Introduction

Erythrocytes or red blood cells (RBCs) are the most numerous type of cells in human blood, and are responsible for the transport of oxygen and carbon dioxide. Typically, for healthy human beings, at a standard altitude, females average about 4.8 million of these cells per cubic millimetre of blood, while males average about 5.4 million per cubic millimetre. The lifespan of RBCs is approximately 120 days. Thereafter, they are ingested by phagocytic cells in the liver and spleen (approx. three million RBCs die and are scavenged each second), and the iron in their hemoglobin (which gives them their characteristic dark colour) is reclaimed for reuse. The remainder of the heme portion of the molecule is degraded into bile pigments and excreted by the liver. The typical bi-concaval shape of RBC is the optimal combination of surface area to volume ratio. This shape also provides unique deformability characteristics to the cell, giving it advantageous properties in order to perform its function in small capillaries. Deviation from the usual healthy cell morphology can lead to a loss of normal function and reduced RBC survival. Hence, measurement of RBC shape is an important parameter to describe RBC function.

A significant part of determining the characteristics of RBC is achieved via optical measurements. Ideally, one would like to perform numerical simulations, in order to minimize time-consuming laboratory tests. Accordingly, the objective of this work is to develop a simple approach to ascertain the light scattering response of large numbers of randomly distributed and oriented red blood cells. Due to the fact that the diameter of a typical red blood cell is on the order of 8 μm (*d*≈8×10^{−6} m), which is much larger than the wavelengths of visible light (3.8×10^{−7} m≤*λ*≤7.2×10^{−7} m), geometric ray-tracing can be used to determine the amount of propagating optical energy, characterized by the Poynting vector, that is reflected and absorbed by multiple RBCs.1 Ray-tracing is highly amenable to rapid large-scale computation needed to track the scattering of incident light beams, comprised of multiple rays, by multiple cells (figure 1), thus making it an ideal simulation paradigm.

The specific model problem that we consider is an initially coherent beam (figure 1), composed of multiple co-linear rays, where each ray is a vector in the direction of the flow of electromagnetic (optical) energy, which, in isotropic media, corresponds to the normal to the wave front. Thus, for isotropic media, the rays are parallel to the wave's propagation vector (figure 1). Of particular interest is to describe the break-up of initially highly directional coherent beams, for example lasers, which do not spread out into multidirectional rays unless they encounter multiple scatterers. The overall objective of the work is to provide a straightforward approach that can be implemented by researchers in the field, using standard desktop computers.

## 2. Parametrization of cell configurations

One of the most widely cited bi-concaval representations for RBC's (figure 1) is (Evans & Fung 1972)(2.1)The outward surface normals, ** n**, needed later during the scattering calculations, are easy to characterize by computing

**=∇**

*n**F*/‖∇

*F*‖. The orientation of the cells, usually random, can be controlled, via standard rotational coordinate transformations, with random angles (figure 1).

The classical random sequential addition algorithm (Widom 1966) is used to place non-overlapping cells randomly into the domain of interest. This algorithm is adequate for the volume fraction range of interest. However, if higher volume fractions are desired, more sophisticated algorithms, such as the equilibrium-based Metropolis algorithm can be used. See Torquato (2002) for a detailed review of such methods. Furthermore, for much higher volume fractions, effectively packing (and ‘jamming’) particles to theoretical limits, a new novel class of methods, based on simultaneous particle flow and growth, has been developed by Torquato and coworkers (see, for example, Kansaal *et al*. (2002) and Donev *et al*. (2004, 2005*a*,*b*)).

## 3. Plane harmonic electromagnetic waves

Recall that, in free space, the propagation of light can be described via an electromagnetic formalism, Maxwell's equations, presented here in simplified form(3.1)and(3.2)where ** E** is the electric field intensity,

**is the magnetic flux intensity,**

*H**ϵ*

_{0}is the permittivity and where

*μ*

_{0}is the permeability. Using standard vector identities, one can show that(3.3)and that(3.4)where the speed of light is . Now consider the case of plane harmonic waves, for example of the form(3.5)where

**is an initial position vector to the wave front and the wavenumber**

*r***is their direction of propagation. For plane waves,**

*k***.**

*k***=constant. We refer to the phase as**

*r**ϕ*=

**.**

*k***−**

*r**ωt*, and

*ω*=2

*π*/

*τ*as the angular frequency, where

*τ*is the period. For ‘plane waves’, the wave front is a plane on which

*ϕ*is constant, which is orthogonal to the direction of propagation, characterized by

**. In the case of harmonic waves, we have(3.6)and**

*k***.**

*k***=0 and**

*E***.**

*k***=0. The three vectors,**

*H***,**

*k***and**

*E***constitute a mutually orthogonal triad. The direction of ray propagation is given by**

*H***×**

*E***/‖**

*H***×**

*E***‖. Since the free space propagation velocity is given by for an electromagnetic wave in a vacuum and for electromagnetic waves in another medium, we can define the index of refraction as2(3.7)**

*H*### 3.1 Optical energy propagation

Light waves travelling through space carry electromagnetic (optical) energy, which flows in the direction of wave propagation. The energy per unit area per unit time flowing perpendicularly into a surface in free space is given by the Poynting vector ** S**=

**×**

*E***. Since at optical frequencies**

*H***,**

*E***and**

*H***oscillate rapidly, it is impractical to measure instantaneous values of**

*S***directly. Now consider the harmonic representations in equation (3.5) which leads to(3.8)and consequently the average value over a longer time interval than the time scale of rapid random oscillation,(3.9)We define the irradiance as(3.10)Thus, the rate of flow of energy is proportional to the square of the amplitude of the electric field. Furthermore, in isotropic media, which we consider for the duration of the work, the direction of energy is in the direction of**

*S***and in the same direction as**

*S***. Since**

*k**I*is the energy per unit area per unit time, if we multiply by the ‘cross-sectional’ area of the ray (

*a*

_{r}), we obtain the energy associated with the ray, denoted as

*Ia*

_{r}=

*Ia*

_{b}/

*N*

_{r}, where

*a*

_{b}is the cross-sectional area of a beam (comprising all of the rays) and

*N*

_{r}is the number of rays in the beam (figure 1).

### 3.2 Reflection and absorption of energy

One appeal of geometrical optics is that elementary concepts are employed. For example, the law of reflection describes how light is reflected from smooth surfaces (figure 2). The angle between the point of contact of a ray and the outward normal to the surface at that point is the angle of incidence (*θ*_{i}). The law of reflection states that the angle at which the light is reflected is the same as the angle of incidence and that the incoming (incident) and outgoing (reflected, *θ*_{r}) rays lay in the same plane and *θ*_{i}=*θ*_{r}. The *law of refraction* states that, if the ray passes from one medium into a second one (with a different index of refraction), and, if the index of refraction of the second medium is less than that of the first, the angle the ray makes with the normal to the interface is always less than the angle of incidence, and can be can be written as (the law of refraction)(3.11)where *θ*_{t} is the angle of the transmitted ray (figure 2).

We consider a plane harmonic wave incident upon a plane tangent to the boundary separating two optically different materials, which produces a reflected wave and a transmitted (refracted) wave (figure 2). The amount of incident electromagnetic energy (*I*_{i}) that is reflected (*I*_{r}) is given by the total reflectance(3.12)where 0≤*R*≤1 and where, for unpolarized (natural) light (see appendix A),(3.13)where is the ratio of the refractive indices of the ambient (incident) medium (*n*_{i}) and transmitted cell medium (*n*_{t}), , where is the ratio of the magnetic permeabilities of the surrounding incident medium (*μ*_{i}) and transmitted cell medium (*μ*_{t}), .

For most materials, the magnetic permeability is, within experimental measurements, virtually the same.3 For the remainder of the work, we shall take , i.e. *μ*_{0}=*μ*_{i}≈*μ*_{t}.

Henceforth, we assume that the medium surrounding the cells behaves as a vacuum, thus, there are no energetic losses as the electromagnetic rays pass through it. Furthermore, we assume that all electromagnetic energy that is absorbed by a cell becomes trapped, and is not re-emitted. This assumption is discussed further later.

## 4. Computational algorithm

The primary quantity of interest is the behaviour of the propagation of the optical energy, characterized by the irradiance. For example, consider the following metrics for overall irradiance of the beam:(4.1)where *N*_{r} is the number of rays comprising the beam and where *I*_{0}=‖*I*(0)‖ is the magnitude of the initial irradiance at time *t*=0. The computational algorithm is as follows, starting at *t*=0 and ending at *t*=*T*:(4.2)

The time step size Δ*t* is dictated by the size of the cells. A somewhat *ad-hoc* approach is to scale the time step size according to Δ*t*∝*ξb*/‖** v**‖, where

*b*is the radius of the cells, ‖

**‖ is the magnitude of the velocity of the rays and**

*v**ξ*is a scaling factor, typically 0.05≤

*ξ*≤0.1.

For step (1), it is convenient to determine whether a ray has just entered a cell domain by checking if , where are the coordinates of the cell expressed in a rotated frame that is aligned with the axes of symmetry of the cell and then to compute the normal ** n**=∇

*F*/‖∇

*F*‖ in that frame.

## 5. A computational example

### 5.1 System parameters

We considered groups of randomly dispersed equal-sized cells, increasing in number, *N*_{c}=1000, 2000, 4000 and 8000, in a rectangular domain of dimensions (figure 3), 1 mm×1 mm×1 cm. This corresponds to a section of a standard testing device, described in detail in §5.2. The stated number of cells corresponded to standard testing hematocrit values. The cells' major diameter was the nominal value of *d*=8×10^{−6} m. A commonly used set of geometric parameters for the cell in equation (2.1) are given by Evans & Fung (1972) as *c*_{0}=0.207161, *c*_{1}=2.002558 and *c*_{2}=−1.122762. The beam was of circular cross-section with diameter of 0.79375 mm (1/32 of an inch, which falls in the range of beams used in experiments described later). The irradiance (Poynting vector magnitude) beam parameter was set to *I*=*I*_{0}*N*−*m*/(*m*^{2}−sec), where the irradiance for each ray was calculated as *I*_{0}*a*_{b}/*N*_{r}, where *a*_{b} was the cross-sectional area of the beam.4 We used successively higher ray densities, *N*_{r}=200, 400, 600, 800 and 1000, etc. rays (figure 3), to represent the beam. The simulations were run until the rays completely exited the domain, which corresponded to a time-scale on the order of (10^{−2} m)/*c*, where *c* is the speed of light. The initial velocity vector for all of the initially co-linear rays comprising the beam was ** v**=(

*c*,0,0).

### 5.2 Computational results

The ratio of the refractive indices varies around 1.0. The exact value corresponds to the state of the cell, including membrane characteristics and haemoglobin concentration. We chose a ratio of refractive indices of , which is consistent with values commonly found in the literature. As the plots in figure 4 indicates, the total amount of energy that is forwardly scattered (defined as the components Poynting ray vectors in the positive *x*-direction) for decreases with the number of cells (scatterers).5 A sequence of frames of the typical ray motion is provided in figure 3. Table 1 tabulates the transmitted energy for various numbers of cells present. *It is important to emphasize that these calculations were performed within a few minutes on a single standard (DELL Precision) laptop*.

Computational tests with higher ray resolution were also performed. We increased the ray density up to 10 000 rays (starting from 200 rays), but found negligible change with respect to the 1000 ray resolution simulation. Thus, beyond *N*_{r}=1000 rays, the computational results changed negligibly, and can be considered to have converged. This cell/ray system provided stable results, i.e. increasing the number of rays and/or the number of cells surrounding the beam resulted in negligibly different overall system responses. Of course, there can be cases where much higher resolution may be absolutely necessary. Thus, it is important to note that a straightforward, natural, algorithmic parallelism is possible with this computational technique. This can be achieved in two possible ways: (i) by assigning each processor its share of the rays, and checking which cells make contact with those rays or (ii) by assigning each processor its share of cells, and checking which rays make contact with those cells.

### 5.3 Laboratory experiments

#### 5.3.1 Preparation of human and murine erythrocytes (RBC)

Blood samples from healthy donors were collected in EDTA anticoagulant, after informed consent, at the Children's Hospital Oakland Research Institute. Whole blood was kept at 4 °C and used within 24 h. RBC were isolated by centrifugation, washed three times in HEPES-buffered saline and the buffy coat was removed after each wash. RBC were re-suspended at 30% haematocrit in HEPES buffered saline (150 mM NaCl, 10 mM HEPES, pH 7.4) and stored at 4 °C until used within 48 h. Before use, cells were suspended in buffer at room temperature to a cell concentration as indicated. The exact cell count in the suspension was determined using the Guava Easycount flowcytometer (Guava Technologies, Hayward, CA).

#### 5.3.2 Light scatter measurements

1.5 ml of cell suspension containing the indicated cell concentration in a cuvet with a 1 cm light path was put in a Varian 50 Cary Bio spectrophotometer (Varian Analytical Instruments, Palo Alto, CA). Light transmittance (*T*=*I*_{x}/‖*I*(0)‖), defined as the ratio of intensity of detected light (*I*_{x}) relative to incoming light (‖*I*(0)‖) of cell suspensions relative to buffer without cells was recorded and averaged over a 1 min interval. Wavelengths were varied from 200 to 800 nm as indicated and specific measurements were performed at 420 and 710 nm, wavelengths of maximum and minimum light absorbance, respectively. In addition, the intensity of the incoming beam was restricted to approximately 1% of the original intensity by a neutral filter.

### 5.4 Comparison between computational predictions and experimental results

In the range of cell concentrations tested, the computational predictions and laboratory results are in close agreement, as indicated in figure 5 and tables 1–3. Although the computations corresponded closely to both wavelengths of light, the match is closer to the 710 nm wavelength, since that wavelength reflects in a manner more consistent with the ratio of refractive indices used in the computations, as opposed to the 420 nm wavelength light which is nearly a purely absorbing combination with RBC.

Figure 5 shows the relative light transmittance *T* as a function of the number of cells per millilitre for different wavelengths of light. Whereas the incoming light (*I*(0)) was greatly affected by placing masks with different circular cross-sections in the light path, the transmittance *T* was not affected. The diameter of 1/32 of an inch for the diameter of the beam used for computation falls within the size used in our experimental approach. Furthermore reducing the incoming light to 1% of it is original value by the use of a neutral filter did not affect the transmittance. The data indicated in figures and tables were collected without restriction on the incoming light. Together, these data indicate that the beam intensity chosen for the computational model corresponded to the experimental approach.

We remark that, in the computations, the refracted energy absorbed by the cells was assumed to remain trapped within the cell. Certainly, some of the absorbed energy by the cells is converted into heat. An analysis of the thermal conversion process can be found in Zohdi (2006). Another level of complexity involves dispersion when light is transmitted through cells. Dispersion is the decomposition of light into its component wavelengths (or colours), which occurs because the index of refraction of a transparent medium is greater for light of shorter wavelengths. Accounting for dispersive effects is quite complex since it leads to a dramatic growth in the number of rays.

## 6. Extensions and concluding remarks

In summary, the objective of this work was to develop a simple computational framework, based on geometrical optics methods, to rapidly determine the light scattering response of multiple red blood cells. Because the wavelength of light (3.8×10^{−7}≤*λ*≤7.2×10^{−7} m) is approximately an order of magnitude smaller than the typical red blood cell scatterer (*d*≈8×10^{−6} m), geometric ray-tracing theory is applicable, and can be used to rapidly ascertain the amount of propagating optical energy, characterized by the Poynting vector, that is reflected and absorbed by multiple cells. Three-dimensional examples were given to illustrate the technique, and the computational results match closely with experiments performed on blood samples at the red cell laboratory in CHORI.

We conclude by stressing a few points for possible extensions. First, a more general way to characterize a wider variety of RBC states, which are not necessarily always bi-concaval, can be achieved by modifying the equation for a generalized ‘hyper’-ellipsoid:(6.1)where the *s*'s are exponents. Values of *s*<1 produce non-convex shapes, while *s*>2 values produce ‘block-like’ shapes. Furthermore, we can introduce particulate aspect ratio, defined by , where *r*_{2}=*r*_{3}, AR>1 for prolate geometries and AR<1 for oblate shapes. To produce the shape of a typical RBC, we introduce an extra term in the denominator of first axes term:(6.2)where and *c*_{1}≥0 and *c*_{2}≥0. The effect of the term, is to make the effective radius of the ellipsoid in the *x*-direction grow as one moves away from the origin. As before, the outward surface normals, ** n**, needed during the scattering calculations, are easy to characterize by writing

**=∇**

*n**F*/‖∇

*F*‖ with respect to a rotated frame that is aligned with the axes of symmetry of the generalized cell.

Second, it is important to recognize that one can describe the aggregate ray behaviour in a more detailed manner via higher moment distributions of the individual ray-fronts and their velocities. For example, consider any quantity, *Q*, with a distribution of values (*Q*_{i}, *i*=1, 2, …*N*_{r}=rays) about an arbitrary reference value, denoted *Q*^{☆}, as follows, , where . The various moments characterize the distribution, for example: (i) measures the first deviation from the average, which equals zero; (ii) is the average; (iii) is the standard deviation; (iv) is the skewness; and (v) is the kurtosis. The higher moments, such as the skewness, measure the bias or asymmetry of the distribution of data, while the kurtosis measures the degree of peakedness of the distribution of data around the average.

Finally, when more microstructural features are considered, for example, topological and thermal variables, parameter studies become quite involved. In order to eliminate a trial and error approach to determine the characteristics of the types of cells that would be needed to achieve a certain level of scattering, in Zohdi (in press), an automated computational inverse solution technique has recently been developed to ascertain scatterer combinations which deliver prespecifed electromagnetic scattering, thermal responses and radiative (infrared) emission, employing genetic algorithms in combination with implicit staggering solution schemes, based upon approaches found in Zohdi (2002, 2003*a*,*b*, 2004*a*,*b*,*c*).

Generally, RBC behaviour under fluid shear stress and response to osmolality changes is essential for normal function and survival. The ability to predict and measure shape and deformation of individual RBC under fluid shear stress will improve diagnosis of RBC disorders and impact new avenues to treatment. New nano-technology approaches coupled with real time computational analysis will make it feasible to generate shape and deformability histograms in very small volumes of blood. This line of research is currently being pursued by the authors, in particular to help detect blood disorders, which are characterized by the deviation of the shape of cells from those of healthy ones under standard test conditions. Such disorders, in theory, could be detected by differences in their scattering responses from that of healthy cells.

Red cell shape is essential for proper function in the circulation. Changes in shape will lead to a decreased red cell survival often accompanied by anaemia. Genetic disorders of cytoskeletal proteins will lead to red cell pathology including hereditary spherocytosis and hereditary eliptocyosis (Eber & Lux (2004) and Gallagher (2004*a*,*b*)). Changes in membrane and cytosolic proteins may affect the state of hydration of the cell and thereby its morphology. Millions of humans are affected by haemoglobinopathies such as Sickle cell disease and Thalassemia (Forget & Cohen (2005) and Steinberg *et al*. (2005)). The altered hemoglobin in these disorders can lead to changes in red cell properties, including membrane damage. Any of these conditions will result in an alteration of the scattering properties of the population of red cells. It is hoped that simple scatter measurements and fitting of the obtained data to our simulation model will reveal altered parameters of the red cell population related to red cell pathology. We hypothesize that this approach may be used as part of the diagnostic process or to evaluate treatment. Changes in clinical care may show a trend to normalization of red cell scatter characteristics, and therefore an improvement of red cell properties.

## Footnotes

↵See Bohren & Huffman (1998), Elmore & Heald (1985) and van de Hulst (1981).

↵All electromagnetic radiation travels at the speed of light in a vacuum,

*c*≈3×10^{8}m s^{−1}. A more precise value, given by the National Bureau of Standards, is*c*≈2.997924562×10^{8}±1.1 m s^{−1}. For visible light, 3.8×10^{−7}m≤*λ*≤7.2×10^{−7}m.↵A few notable exceptions are concentrated magnetite, pyrrhotite, and titanomagnetite (Telford

*et al*. (1990) and Nye (1957)).↵Because of the normalized structure of the metric, it is insensitive to the magnitude of

*I*_{o}for the scattering calculations. The initial magnitude of the Poynting vector magnitude is , where, initially, only one component is nonzero,*I*_{x}(0)=*I*_{o}, in the*x*direction.↵The system at time

*t*=*T*indicated that all rays have exited the scattering system.↵Throughout the analysis we assume that .

↵The limiting case , is the critical angle (

*θ*_{i}^{*}) case.- Received April 22, 2006.
- Accepted June 1, 2006.

- © 2006 The Royal Society