Discrimination of granulocyte subtypes from light scattering : theoretical analysis using a granulated sphere model

We perform extensive simulations of light scattering by a granulated sphere in the size and refractive index range of human granulated leucocytes using the discrete dipole approximation. We calculate total and depolarized side scattering signals as a function of the size and refractive indices of cell and granules, and the granule volume fraction. Using typical parameters derived from the literature data on granulocyte morphology, we show that differences between experimentally measured signals of two granulocyte subtypes can be explained solely by the difference in their granule sizes. Moreover, the calculated depolarization ratio quantitatively agrees with experimental results. We also use the Rayleigh-Debye-Gans approximation and its second order extension to derive analytical expressions for side scattering signals. These expressions qualitatively describe the scaling of signals with varying model parameters obtained by rigorous simulations, and even lead to quantitative agreement in some cases. Finally, we show and discuss the dependence of extinction efficiency and asymmetry parameter on size and volume fraction of granules. ©2007 Optical Society of America OCIS codes: (000.4430) Numerical approximation and analysis; (170.1530) Cell analysis; (290.4210) Multiple scattering; (290.5850) Scattering, particles; (290.5855) Scattering, polarization. References and links 1. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles, (Cambridge University Press, Cambridge, 2002). 2. P. Chylek, G. Videen, D. J. W. Geldart, J. S. Dobbie, and H. C. W. Tso, "Effective medium approximations for heterogeneous particles," in Light Scattering by Nonspherical Particles, Theory, Measurements, and Applications, M. I. Mishchenko, J. W. Hovenier, and L. D. Travis, eds. (Academic Press, New York, 2000), pp. 273-308. 3. L. Kolokolova and B. A. S. Gustafson, "Scattering by inhomogeneous particles: microwave analog experiments and comparison to effective medium theories," J. Quantum Spectrosc. Radiat. Transfer 70, 611625 (2001). 4. N. V. Voshchinnikov, V. B. Il'in, and T. Henning, "Modelling the optical properties of composite and porous interstellar grains," Astron. Astrophys. 429, 371-381 (2005). 5. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Multiple Scattering of Light by Particles: Radiative Transfer and Coherent Backscattering, (Cambridge University Press, Cambridge, 2006). 6. Y. L. Xu and B. A. S. Gustafson, "Comparison between multisphere light-scattering calculations: Rigorous solution and discrete-dipole approximation," Astrophys. J. 513, 894-909 (1999). 7. D. W. Mackowski and M. I. Mishchenko, "Calculation of the T matrix and the scattering matrix for ensembles of spheres," J. Opt. Soc. Am. A 13, 2266-2278 (1996). 8. M. I. Mishchenko, L. Liu, D. W. Mackowski, B. Cairns, and G. Videen, "Multiple scattering by random particulate media: exact 3D results," Opt. Express 15, 2822-2836 (2007). 9. A. Taflove and S. C. Hagness, Advances in Computational Electrodynamics: the Finite-Difference TimeDomain Method, 3rd ed., (Artech House, Boston, 2005). #88472 $15.00 USD Received 10 Oct 2007; revised 12 Nov 2007; accepted 12 Nov 2007; published 29 Nov 2007 (C) 2007 OSA 10 December 2007 / Vol. 15, No. 25 / OPTICS EXPRESS 16561 10. M. A. Yurkin and A. G. Hoekstra, "The discrete dipole approximation: an overview and recent developments," J. Quantum. Spectrosc. Radiat. Transfer 106, 558-589 (2007). 11. K. Lumme and J. Rahola, "Light-scattering by porous dust particles in the discrete-dipole approximation," Astrophys. J. 425, 653-667 (1994). 12. A. K. Dunn, "Modelling of light scattering from inhomogeneous biological cells," in Optics of Biological Particles, A. G. Hoekstra, V. P. Maltsev, and G. Videen, eds. (Springer, London, 2006), pp. 19-29. 13. Flow Cytometry and Sorting, 2nd ed., M. R. Melamed, T. Lindmo, and M. L. Mendelson, eds. (Wiley-Liss, New York, 1990). 14. B. G. de Grooth, L. W. Terstappen, G. J. Puppels, and J. Greve, "Light-scattering polarization measurements as a new parameter in flow cytometry," Cytometry 8, 539-544 (1987). 15. S. Suzuki and N. Eguchi, "Leukocyte differential analysis in multiple laboratory species by a laser multiangle polarized light scattering separation method," Exp. Anim. 48, 107-114 (1999). 16. M. Hedhammar, M. Stenvall, R. Lonneborg, O. Nord, O. Sjolin, H. Brismar, M. Uhlen, J. Ottosson, and S. Hober, "A novel flow cytometry-based method for analysis of expression levels in Escherichia coli, giving information about precipitated and soluble protein," J. Biotech. 119, 133-146 (2005). 17. F. Lavergne-Mazeau, A. Maftah, Y. Cenatiempo, and R. Julien, "Linear correlation between bacterial overexpression of recombinant peptides and cell light scatter," Appl. Environ. Microbiol. 62, 3042-3046 (1996). 18. L. W. M. M. Terstappen, B. G. de Grooth, K. Visscher, F. A. van Kouterik, and J. Greve, "Four-parameter white blood cell differential counting based on light scattering measurements," Cytometry 9, 39-43 (1988). 19. S. Lavigne, M. Bosse, L. P. Boulet, and M. Laviolette, "Identification and analysis of eosinophils by flow cytometry using the depolarized side scatter-saponin method," Cytometry 29, 197-203 (1997). 20. S. L. Perkins, "Normal blood and bone marrow values in humans," in Wintrobe's Clinical Hematology, 11th ed., J. P. Greer, J. Foerster, and J. N. Lukens, eds., (Lippincott Williams & Wilkins Publishers, Baltimore, USA, 2003), pp. 2738-2741. 21. K. M. Skubitz, "Neutrophilic leukocytes," in Wintrobe's Clinical Hematology, 11th ed., J. P. Greer, J. Foerster, and J. N. Lukens, eds. (Lippincott Williams & Wilkins Publishers, Baltimore, USA, 2003), pp. 267-310. 22. P. Lacy, A. B. Becker, and R. Moqbel, "The human eosinophil," in Wintrobe's Clinical Hematology, 11th ed., J. P. Greer, J. Foerster, and J. N. Lukens, eds., (Lippincott Williams & Wilkins Publishers, Baltimore, USA, 2003), pp. 311-334. 23. H. P. Ting-Beall, D. Needham, and R. M. Hochmuth, "Volume and osmotic properties of human neutrophils," Blood 81, 2774-2780 (1993). 24. P. Brederoo, J. van der Meulen, and A. M. Mommaas-Kienhuis, "Development of the granule population in neutrophil granulocytes from human bone marrow," Cell Tissue Res. 234, 469-496 (1983). 25. L. W. Diggs, D. Sturm, and A. Bell, The Morphology of Human Blood Cells, 5th ed., (Abbott Laboratories, Abbott Park, IL 60064, 1985). 26. G. J. Puppels, H. S. P. Garritsen, G. M. J. Segers-Nolten, F. F. M. de Mul, and J. Greve, "Raman microspectroscopic approach to the study of human granulocytes," Biophys. J. 60, 1046-1056 (1991). 27. S. A. Livesey, E. S. Buescher, G. L. Krannig, D. S. Harrison, J. G. Linner, and R. Chiovetti, "Human neutrophil granule heterogeneity: immunolocalization studies using cryofixed, dried and embedded specimens," Scanning Microsc. Suppl 3, 231-239 (1989). 28. W. T. Daems, "On the fine structure of human neutrophilic leukocyte granules," J. Ultrastruct. Res. 24, 343348 (1968). 29. O. W. Bjerrum, "Human neutrophil structure and function with special reference to cytochrome b559 and beta 2-microglobulin," Dan. Med. Bull. 40, 163-189 (1993). 30. D. F. Bainton, "Neutrophilic leukocyte granules: from structure to function," Adv. Exp. Med. Biol. 336, 1733 (1993). 31. A. E. Zharinov, P. A. Tarasov, A. N. Shvalov, K. A. Semyanov, D. R. van Bockstaele, and V. P. Maltsev, "A study of light scattering of mononuclear blood cells with scanning flow cytometry," J. Quantum. Spectrosc. Radiat. Transfer 102, 121-128 (2006). 32. K. A. Semyanov, A. E. Zharinov, P. A. Tarasov, M. A. Yurkin, I. G. Skribunov, D. R. van Bockstaele, and V. P. Maltsev, "Optics of leucocytes," in Optics of Biological Particles, A. G. Hoekstra, V. P. Maltsev, and G. Videen, eds. (Springer, London, 2006), pp. 253-264. 33. C. F. Bohren and D. R. Huffman, Absorption and scattering of Light by Small Particles, (Wiley, New York, 1983). 34. 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. Transfer 106, 546-557 (2007). 35. "Amsterdam DDA," http://www.science.uva.nl/research/scs/Software/adda (2007). 36. "Description of the national compute cluster Lisa," http://www.sara.nl/userinfo/lisa/description/ (2005). 37. N. W. Ashcroft and J. Lekner, "Structure and resistivity of liquid metals," Phys. Rev. 145, 83-90 (1966). 38. M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids, (Oxford University Press, Oxford, 1989). 39. V. V. Tuchin, L. V. Wang, and D. A. Zimnyakov, Optical Polarization in Biomedical Applications, (Springer, Berlin, 2006). #88472 $15.00 USD Received 10 Oct 2007; revised 12 Nov 2007; accepted 12 Nov 2007; published 29 Nov 2007 (C) 2007 OSA 10 December 2007 / Vol. 15, No. 25 / OPTICS EXPRESS 16562


Introduction
Simulation of light scattering by homogeneous particles with simple shapes, such as spheroids, is relatively easy [1]. Unfortunately, most natural particles are inhomogeneous and have a complex shape. A broad class of such particles can be characterized by a matrix with multiple inclusions (granules). A number of approximations exist to simulate light scattering by these particles which fall in one of two connected realms: effective medium theories (EMTs) [2][3][4] or radiative transfer theories [5]. The latter is especially relevant when an infinite medium is considered or granule positions randomly fluctuate with significant amplitude during the observation time.
However, in many cases none of these approximations is accurate enough and rigorous light scattering theories are required. The superposition T-matrix method allows simulation of light scattering by clusters of spheres [6,7], which is a particular case of a granulated particle when the matrix is vacuum. Recently Mishchenko et al. [8] studied clusters of spheres with respect to radiative transfer theories. Apart from that, only general methods like the finite difference time domain method (FDTD) [9] or the discrete dipole approximation (DDA) [10] are capable of simulating light scattering by arbitrarily shaped granulated particles. Most researchers who studied granulated particles are concerned with astrophysical or atmospheric applications [3,4,8,11], while our main application -light scattering by biological cells -is a much less studied field.
Biological cells, when suspended in liquid, have an important advantage with respect to the light scattering simulation. Their relative refractive index is close to unity. This accelerates the rigorous methods and improves the accuracy of the EMTs and other approximate theories. FDTD simulation of light scattering by biological cells was performed in a number of manuscripts by Dunn and coworkers (summarized in [12]). To the best of our knowledge, these are the only studies that include cell organelles and granules in rigorous light scattering modeling. The main conclusion is that side scattering by biological cells are mainly determined by small granules, which agrees with the general notion of flow cytometry [13]. However, only the dependence on the volume fraction and not on the size of the granules was studied.
Study of light scattering by granulated biological cells is motivated by several empirical techniques that use the dependence of light scattering signals measured by flow cytometer on the characteristics of the granules [14][15][16][17]. All these techniques are successfully used in practice but lack a rigorous theoretical foundation. For instance, production of the inclusion bodies consisting of protein synthesized by E. coli cells was found to correlate with forward scattering [16,17], which can be qualitatively explained by the increase of total amount of matter that scatters light inside the cells. Probably, the most fascinating application of light scattering to characterize granulated cells are the results of de Grooth, Terstappen and coworkers [14,18,19], who showed that neutrophils and eosinophils can be discriminated by their depolarized side scattering intensity I ⊥ . This result was extended to other species by Suzuki and Eguchi [15]. The values for light scattering signals are given in arbitrary units by de Grooth et al. [14] and, hence, can be compared to other studies only qualitatively. However, the ratios of the signals, which are also provided, can be compared quantitatively to predictions of any models, thus providing a stringent test of the latter.
Granulocytes are the most numerous type of leukocytes, consisting of three subtypes: neutrophils, eosinophils, and basophils. Their abundance relative to all leukocytes in normal condition is in the ranges 46-73%, 0-4.4%, and 0.2-1.2% respectively [20]. The major role of neutrophils is to protect the host against infectious agents. To accomplish this task, the neutrophil must first sense infection, migrate to the site of the infecting organism, and then destroy the infectious agents [21]. Although rare in healthy individuals, the eosinophil is prominent in peripheral blood and tissue in association with various disease conditions including allergy, inflammatory responses against metazoan helminthic parasites, and certain skin and malignant conditions [22].
In this manuscript we aim at two goals: to get general insight into light scattering by index-matching granulated particles much larger than the wavelength and to rigorously explain why eosinophils and neutrophils have different I ⊥ . For that we construct a simple model of granulocytes in Section 2, starting from literature data on its morphology. In Section 3 we specify the total and depolarized side scattering signals measured by flow cytometers. We derive analytical expressions for these signals using the Rayleigh-Debye-Gans (RDG) approximation and its second-order extension in Sections 4 and 5 respectively. We show the results of extensive DDA simulations in Section 6, comparing them with predictions of approximate theories and known experimental results. We conclude the manuscript in Section 7 and discuss perspectives for characterizing granulated particles, specifically granulocytes, with light scattering.

Simple granulocyte model and DDA simulations
Both neutrophils and eosinophils have a size of 10-15 μm as measured in blood smears [21]. However Ting-Beall et al. [23] reported that a normal suspended neutrophil has a volume of 299 ± 64 μm 3 , which corresponds to a diameter 8.3 ± 0.6 μm (intervals correspond to two standard deviations). A similar value for the diameter of mature neutrophils (8.3 ± 1.2 μm) was obtained by Brederoo et al.
[24] using electron microscopy. For eosinophils a value of 8 μm has been reported [22]. The nucleus is separated into definite lobes with a very narrow filament or strand connecting the lobes. Usually, two or three lobes are observed [22,25,26].
Neutrophils contain two major types of granules: azurophilic and specific granules, although more minor classes can be specified [24,27,28]. Azurophilic granules are larger and more dense than the specific ones, and the ratio of their numbers is 1/2-1/3 [29]. A typical size of azurophilic granules is 500 nm, whereas specific granules are either spherical with a size of 200 nm (Livesey, et al., [27] reported the range 100-250 nm) or rod shaped (130×1000 nm) [30]. On average a cell cross section contains 200-300 granules [30], the total number per whole cell is 1900-6300, with a total volume of 30 μm 3 [29]. Eosinophils contain three major classes of granules: crystalloid, primary, and small granules. Crystalloid granules measure 0.5-0.8 μm in diameter [22] (Puppels et al. [26] reported another range: 0.9-1.3 μm) and contain crystalline electron-dense cores surrounded by an electron-lucent matrix. There are approximately 200 crystalloid granules in each cell [22,26]. Primary granules measure 0.1-0.5 μm in diameter and are less abundant than crystalloid granules. Small granules are even smaller than primary ones [22].
We use a granulated sphere as a model for granulocytes (see Fig. 1). The sphere contains no nucleus and all the granules are identical and randomly placed inside the cell. Thus, we compromise between realistic shapes (cf. [12]) and a general fundamental approach (cf. [8]). We are able to study the dependence of light scattering of each of the particle parameters independently, and at the same time to quantitatively describe at least the depolarization ratios of neutrophils and eosinophils (see Section 6). In the following we describe the default parameters of the model. The diameter of the cell D c = 8 μm and volume fraction of granules f = 0.1 correspond to the above data on morphology of both neutrophils and eosinophils. Dunn collected information on refractive indices of the cell cytoplasm and constituents from several sources [12]. However, the ranges are broad, which makes it hard to select a particular value. We set the cytoplasm refractive indices, relative to the medium (for which we consider buffered saline m 0 = 1.337), to m c = 1.015, which corresponds to the values obtained by fitting experimental light scattering patterns of lymphocytes by a multi-layered sphere model [31,32]. For the granules we set the relative refractive index to the upper limit as if granules consist entirely of proteins, m g = 1.2. This choice exaggerates all light scattering effects of granules; however, as we will show, all conclusions also are valid for smaller m g . We completely neglect absorption in this study, and set the wavelength to that of a 0.66 μm semiconductor laser, which is 0.4936 μm in the medium. The main variable is the granule diameter d g , which is varied from 0.075 to 2 μm. The lower limit is determined by the size of the dipole in DDA, which was set to 0.041 μm for all simulations, corresponding to 12 dipoles per wavelength in the outer medium. Each DDA simulation depends on random granule placement, therefore we repeat it 10 times and show mean value ± 2×SD (standard deviation) for each simulated value. In addition, we simulated the result for effectively d g = 0, by using the Mie theory applied to a homogenous sphere, where the refractive index is obtained by Maxwell-Garnett EMT [33]. Since all refractive indices are close to unity, the differences between most EMTs are negligible.
We also varied some of the parameters from their default values. To limit the number of DDA simulations and keep the time of this study feasible we varied each parameter with others fixed. For each set of parameters a whole range of d g was calculated. We As a numerical implementation of the DDA we have used the ADDA computer code v.0.76, which is capable of running on a cluster of computers (parallelizing a single DDA computation), allowing simulating light scattering by scatterers much larger than a wavelength [34,35]. In this manuscript we use the default ADDA settings for dipole polarizability (lattice dispersion relation), iterative method (quasi minimal residual method), and convergence threshold (relative residual norm less than 10 −5 ) and employ the built-in granule generator. For each particle we calculated Q ext , <cosθ >, and the whole Mueller matrix for the whole scattering angle with steps of 0.5° and 5° in polar θ and azimuthal ϕ angles respectively. Moreover, we refined the grid to 0.5° step in ϕ near the side scattering direction to accurately compute the signals described in Section 3. All simulations were run on the Dutch compute cluster LISA [36]. Typical simulation time is half an hour on 8 nodes (each node has dual Intel Xeon 3.4 GHz processor with 4 GB RAM), but it is about 1.5 hours on 16 nodes for D c = 14 μm.

Orthogonal light scattering
We simulate the side scattering signals exactly as described by de Grooth et al. [14]. The incident light propagates along the z-axis and is polarized along the x-axis, and the particle is located in the origin. The exact side scattering direction is along the y-axis. A lens focused on the particles projects the scattered light on the photodetector; both the lens and the photodetector are centered around and orthogonal to the y-axis. A rectangular diaphragm placed before the lens limits the collection scattering angles to θ = 90° ± Δθ and ϕ = 90° ± Δϕ.
The default values are Δθ = Δϕ = 25°. The measured signals are the total side scattering intensity collected by the photodetector I SS and the total depolarized intensity I ⊥ , when a polarizer with polarization axis along the z-axis is placed before the photodetector. The depolarization ratio is defined as D SS = I ⊥ /I SS .
In the original paper [14] an expression of I SS and I ⊥ in terms of the amplitude scattering matrix (S i ) is derived, however it has a typographical error in one sign. Therefore, we present the corrected derivation in the following and also derive an expression in terms of the Mueller scattering matrix (S ij ). We denote the radius vector of a point on the lens as r, and scattering direction n = r/r, which corresponds to scattering angles θ and ϕ. We introduce unity vectors e ⊥ and e || orthogonal and parallel to the scattering plane respectively, which are the local basis for scattered wave in point r, i.e. {e ⊥ , e || , n} is a right-handed orthogonal basis (according to Bohren and Huffman [33]). These vectors are connected to the basis vectors of a spherical coordinate system by e || = e θ , e ⊥ = −e ϕ . We assume that incident light has unity amplitude. Then two components of the scattered electric field in point r (before the lens) is given by: The lens rotates the propagation vector of the scattered wave from n to e y together with the electric field vector. We introduce the following auxiliary unity vectors: , , | | ) ( . cos sin , sin cos It is easy to show from geometrical considerations that , cot sec tan , cot sec tan The corresponding angle-resolved total and depolarized intensities, which are sensed by the detector, are . , Using Eqs. (1), (4), and (7) and the definition of the Mueller scattering matrix in terms of the amplitude scattering matrix [33] one can obtain: The same result can be obtained using the Mueller matrix formalism.
The final signals, measured by the detector, are proportional to the intensities integrated over the scattering aperture: ). , where all Mueller matrix elements are considered at exact side scattering direction. The particular value of κ is not relevant since we will always consider either qualitative behavior of plots including I SS and I ⊥ or relative values. Next, we summarize the experimental conclusions of de Grooth, et al. [14], that we will explain by rigorous DDA simulations in Section 6. The scatter plot of the sample containing eosinophils and neutrophils in I ⊥ versus I SS coordinates shows clear discrimination of two subtypes. The I SS signals are almost the same, while I ⊥ are larger for eosinophils. Moreover, there is significant correlation between I ⊥ and I SS for both neutrophils and eosinophils. Discrimination between two subtypes is the most pronounced using D SS , with mean values of 0.044 and 0.013 for eosinophils and neutrophils respectively (averaged over a single sample). They also showed that D SS is almost constant when Δϕ is decreased at fixed Δθ for both granulocyte subtypes. Thus, the depolarization signals of granulocytes are inherent in contrast to the aperture depolarization, e.g. of a sphere, which is caused by a finite value of the detector aperture around the exact side scattering direction and decreases to zero with decreasing Δϕ.

The Rayleigh-Debye-Gans approximation
Consider our model of a granulated sphere (see Section 2). The position of the granule centers relative to the origin, which is placed in the cell center, is given by vectors r i , where i is from 1 to N, and N is total number of granules. V c , x c and V g , x g are the volume and the size parameter of the cell and one granule respectively, and the volume fraction of granules is given by 3 . We assume that granules are randomly positioned inside the cell, and they are not allowed to overlap. Furthermore, to simplify our derivations we assume that x g < < x c . It is important to note that x g can be both small and large compared to unity. We formulate the random granule position as: r i is uniformly random inside the sphere with size parameter x c − x g . Thus we neglect the boundary effects due to the non-overlapping of granules, which are significant only for large f in the layer of width of order of d g from the cell boundary. It is important to note that the pair distribution function of granules, which will be discussed below, is much more sensitive to the non-overlapping condition than the probability distribution function of a single granule position. The incident wave propagates along the zaxis and the scattering direction n is described by angles θ and ϕ. We consider the scattering problem of a particle in vacuum, i.e. we divide all relevant quantities by the refractive index of the host medium.
To understand the dependence of side-scattering intensity on the granule size we consider the RDG approximation [33]. We consider the limit of small refractive indices, i.e. |m − 1| < < 1 for both cell cytoplasm and granules. Actually, RDG is strictly valid only when x|m − 1| < < 1, which, strictly speaking, is not true for our particles. We will see, however, that it is capable of describing some features that are obtained by accurate DDA simulations.
Our goal is to derive simple analytical expression rather than to keep all derivations as accurate as possible. We can derive the final RDG result numerically with any required accuracy for any particular granule configuration and then perform numerical averaging over all possible configurations. This is a laborious task, albeit much faster than the DDA simulations. However, such brute-force numerical simulations of RDG-based theory are anyway not accurate enough. Therefore, we prefer to additionally sacrifice some accuracy to derive analytical expressions, which give added value compared to rigorous DDA simulations. This value consists in physical insight into the light scattering problem, e.g. scaling laws, and the opportunity for an approximate solution of the inverse light scattering problem. Although we discuss all employed assumptions and approximations, the accuracy of the final expressions can only be determined by comparison with the DDA simulations.
According to the RDG theory, only the diagonal elements of the amplitude scattering matrix are nonzero [33]: where the particle is divided into N + 1 domains: i = 0 corresponds to the cytoplasm and the rest to N granules. m i and V i are refractive index and volume of each domain. h(V,n) is a form factor given by where we have introduced q = k(e z − n). The form factor for a sphere, which center is in the origin, can be obtained analytically [33]: where r and x are the radius and the size parameter of a sphere, and there is no dependency on the azimuthal angle. The asymptotic behavior of Eq. (14) is: Using the linearity of Eq. (12) in the factor m − 1, we consider separately the spherical cytoplasm with factor m c − 1 and superimposed granules with factor m g − m c . Then, Eq. (12) can be rewritten as where ξ(N) holds all the dependency on granule positions: Averaging of ξ(N) over all possible granule positions is performed independently for each element of the sum leading to the same integral as the one in Eq. (13), therefore ). , The second moment of the absolute value of ξ(N) is given by: Now we consider a particular case by assuming that r i and r j are independent. In other words, we neglect the effect of the non-overlapping condition on the statistical properties of granule positions. This is only valid for sufficiently small volume fractions ( (18) and (20) one can obtain Equation (21) is obtained under the assumptions x g < < x c and 1 < < f (and the RDG approximation itself). In particular, if the limit x g → 0 and N → ∞ is taken, keeping f constant, then  N) is about an order of magnitude smaller than the second. Therefore, we may completely neglect the cytoplasm except for its total volume and the effect it has on the refractive index of the granules. By that we introduce an assumption x c > > 1, and use it for a new evaluation of Eqs. (18) and (19) discarding the condition 1 < < f . Granules are uniformly distributed inside a large volume with volume fraction f, without overlapping. Taking the limit of infinitely large volume (N → ∞, f = const) and finite θ (θ ∼ 1), that is equivalent to a model of a hard spheres liquid, for which it is known that [37,38] ); .Employing the assumptions discussed above together with Eqs. where Eq. (12) was used. Since averaging and integration can be interchanged, First, we analyze the scaling behavior of g f (v). One may show that and finally: where constants ap 2 , 1 C are determined by the aperture angles Δθ and Δϕ. ap 2 C also weakly depends on x, but this is ignored in this discussion.

The second-order Born approximation
To theoretically approach the depolarized intensity we employ the second order of the RDG or, more precisely, the second-order Born (2-Born) approximation. We use the same definitions as given in Section 4, and the same philosophy of preferring simplicity to the best possible accuracy. The internal field inside the scatterer is a sum of two orders (see e.g. [10]): where E inc (r) and E(r) are the incident and total electric field, χ(r) = (m 2 (r) − 1)/4π is the susceptibility of the medium, V 0 is an infinitesimally small spherical exclusion volume around r. ) , ( r r G ′ is the free space dyadic Green's function, defined as where I is the identity dyadic, k = ω/c is the free space wave vector, R = r − r′, R = |R|, and R Rˆ is a dyadic defined as Most of the scattering quantities can be obtained from the scattering amplitude [33]: To simplify our derivations we consider the depolarized intensity Like in the RDG derivations, we divide the scatterer into the whole volume V c with susceptibility χ c and superimposed granules with total volume V grs and susceptibility χ g − χ c ).
The integral in Eq. (39) is decomposed into four parts: a double integral over V c , a double integral over V grs and two cross-terms (over V c and V grs ). The first term is just a second-order Born approximation for a homogenous sphere, which identically equals zero because of symmetry (S 3 is always zero for a sphere). The cross-terms are very similar and they both can be thought of as RDG scattering by V grs for an incident field produced by the RDG applied to V c . One can show that the RDG applied to an index-matching homogenous sphere much larger than the wavelength produces an internal field parallel to the incident one, except at a distance of order of wavelength from the boundary. The x-component of the internal field does not contribute to the depolarization at exact side scattering direction. Therefore, the relative contribution of the cross-terms in Eq. (39) is at least an order of (m c − 1)/[ fx c (m g − m c )] smaller than the remaining granule-granule interaction, and can be neglected. Moreover, interaction of a granule with itself results in identically zero polarization. Hence, where * denotes complex conjugate, R ij = r j − r i , and integration is performed over a granule positioned in the origin. If all four indices: i, j, i′, and j′ are different then R ij and R i′j′ are independent. We assume that R ij and R ij′ ( j ≠ j′) are also independent; by that we neglect the triple correlations. In case arguments of two Green's tensors in Eq. (41) are independent, they can be averaged independently. The result is then zero because changing sign of the x component of R ij inverts the sign of G zx (R ij ) while not affecting any of the exponents in Eq. (41). Therefore, we need to consider only the terms with either (1) i = i′ and j = j′ or (2) i = j′ and j = i′. All the instantiations of (1) are equivalent, and so are the instantiations of (2 where N > > 1 is assumed. To keep the derivations analytical in calculation of H(R ij ) we assume that R ij > > d g . However, the following derivation is expected to be approximately correct also for small R ij , resulting in the correct scaling laws but with different constants. Employing the assumption of distant granules we obtain where P(R) is a radial distribution function of relative granule positions, defined so that , and h Ω is the following result of averaging over the whole solid angle, since all n are equiprobable, The exponent in Eq. (42) was replaced by a cosine function because the imaginary part of h Ω is zero due to the symmetry properties of the integrand. It is easy to show that P(R) for point granules in a sphere of diameter D c is given by independent of f. Taking into account a finite size of granules, which is assumed much smaller than D c , modifies P(R) only for R close to 0 and D c . Moreover, we neglect the boundary effects for large R, because the contribution to the integral in Eq. (46) from a small interval near D c , which has width of order d g , is relatively small. The correction for small R can be obtained in the limit of infinitely large D c , i.e. considering a granulated sphere as an infinite hard sphere liquid (x c > > 1 and x c > > x g ). The structure factor given by Eq. (24) is basically a Fourier transform of P(R) [37], hence P(R) can be obtained by the inverse Fourier transform This integral cannot be computed in closed analytical form. Therefore, we expand S f (q) in series of f and leave only terms up to the second order, since higher orders anyway cannot be computed accurately without consideration of triple-and higher scattering orders. After that one may obtain the final result for P(R): Using Eqs. (45), (47), and (50), the integral in Eq. (46) can be computed; however, it is still cumbersome. To simplify, we notice that h Ω only weakly depends on R, while different parts of h G (R) exhibit very different scaling with respect to R [cf. Eq. (45)]. The first two terms increase when R → 0 [after multiplication by P(R)], and the integral is determined by a small interval of R, up to several d g . On the contrary, the third part of h G (R) has a smooth behavior, and it's integral is determined by the whole range of R. In particular, when evaluating the integral of the third part of h G (R) the cosine term in h Ω can be neglected, since for kR > > 1, which is valid for most of the R range, the contribution of this term is negligible compared to the unity addend. We denote the resulting h Ω as h Ω (x,∞), and one may show that Next we calculate the contribution of the first two parts of h G (R). As we will show below, this contribution is significant only for small x g . Therefore, we expand h Ω (x g ,R) in series of x g , similar to Eq. (51), assuming kR = O(x g ): It is important to note, that Eq. (55) is not completely accurate, because it is determined by R ∼ d g , for which derivation of Eq. (44) is not accurate. Comparing Eq. (55) to Eqs. (53) and (51) one can see that the contribution of the third part of h G (R) is negligible for very small x g but becomes dominating, compared to other parts, for x g larger than a few times x . This conclusion can be derived rigorously for the whole range of x g , using Eq. (52) instead of Eq. (54) for evaluation of the integral in Eq. (55). Physically speaking, for very small granules the depolarized intensity is determined by the short-range interaction of granules that happen to be close to each other. Starting from the main contribution is from the long-range interaction of all granules. Therefore, we may use Eq. (55) for x g up to some arbitrary chosen constant of order less than unity (we choose it to be unity), and set its contribution to zero for larger x g . Fortunately, this exactly corresponds to the assumptions that were used in obtaining Eq.  )  2  ln  522  373  (  52  35  2  ln  6  1  4  1  ) , Using Eqs. (29) and (56) one may derive the scaling behavior of depolarization ratio: Although our derivations do predict O( f 2 ) corrections of D SS for small x g , this correction is not accurate because of higher-order scattering effects which are of the same order. The factor ) ln is almost constant for a range of large x g that we studied, so we may consider both limiting values of D SS to be independent of x g . Hence, our approximate derivations based on the RDG and the 2-Born predicts the step-wise behavior of D SS (x g ).
The 2-Born also can be used to refine the RDG results for I SS . Although this definitely improves the accuracy, it also makes the final result significantly more complex. Derivation for the exact side scattering direction is more cumbersome than that for I ⊥ , because less terms in intermediate equations cancel out. Averaging over the side scattering aperture also does not seem feasible in closed form. Therefore this refinement contradicts our philosophy for approximate theories and we do not pursue it further.

Results and discussion
The results of DDA simulations of side scattering by granulated spheres with the default set of parameters and varying granule diameter is shown in Fig. 2. For each d g mean values of I SS and I ⊥ and error bars corresponding to two SDs are shown. Labels near some of the points indicate the values of d g . One can see that there are two ranges of d g : from 0.1 to 0.25 μm and from 0.4 to 2 μm, which correspond to distinct regions in Fig. 2. The range from 0.25 to 0.4 μm is intermediate, where I SS and I ⊥ strongly albeit differently depend on d g . The two regions in Fig. 2 qualitatively correspond to the neutrophils and eosinophils on the plots by de Grooth et al. [14]. Considering the morphological characteristics of neutrophil and eosinophils granules, the difference in depolarized side scattering can be explained solely by the difference in d g . However, it is important to note, that direct comparison cannot be performed because the model is a significant simplification of the granulocyte morphology. Apart from neglecting a nucleus, it assumes the same size for all granules, while granule population of a real granulocyte is always heterogeneous both in size and refractive index.
The results for several volume fractions are shown in Fig. 3, separately for I SS , I ⊥ , and their ratios. One may note that EMT result (d g = 0) for D SS seems to be different from the limiting value of the DDA simulations [ Fig. 3(c)]. This is explained by the fact that both I SS and I ⊥ are mostly determined by the granules even for the smallest d g = 0.075 μm that we tried. However, for much smaller granules the signals are those of the homogeneous sphere, which has a different depolarization mechanism (see Fig. 5 and Section 5). Therefore, both I SS and I ⊥ decrease monotonically with d g , while D SS "switches" from the limiting value of small granules to that of a homogeneous sphere, which are naturally different.
Apart from this feature, the general behavior of I SS , I ⊥ , and D SS versus d g , is similar to that predicted by Eqs. (29), (56), and (58): both I SS and I ⊥ rapidly increase with d g for small d g and relatively slowly decay for larger d g . D SS behaves in a step-like manner, switching from a plateau for small d g to another one for larger d g . We include the literature data on average values of D SS for neutrophils and eosinophils [14] in Fig. 3(c). One can see that neutrophils approximately correspond to small d g limit for f = 0.1, while eosinophils have the D SS values between the large and small d g limits for the same f. The latter can be explained by the fact that eosinophils have large (> 0.4 μm), small (< 0.25 μm), and intermediate granules, and f = 0.1 is their typical total volume fraction. It is important to note, that this explanation is Although our derivations do provide such corrections they are not expected to be accurate and therefore are not further discussed (see Sections 4 and 5). One also can see that the slope of the curves in Fig. 4 is almost independent on f for small d g , hence D SS is proportional to f in this regime. This scaling law holds with much better accuracy than that for either I ⊥ or I SS .
We show the results of varying the aperture in Fig. 5. We keep Δθ fixed and vary Δϕ from 0° to 25°. The depolarization ratio does depends on Δϕ over the whole range of d g , however this dependence is uniform in the sense, that variation of Δϕ do not change the general behavior of D SS versus d g . In other words, the aperture part of the depolarization of the granulated sphere is significant, especially for large d g , however it is always less than the inherent part. There is one natural exception to the above conclusion -the EMT result for d g = 0, which scales approximately as Δϕ 2 [14].
A plot of I ⊥ versus I SS for several D c is shown in Fig. 6 come close to each other for large d g . That means that the simple scaling derived from the approximate theories works satisfactory. Moreover, this scaling law is even more accurate for D SS . The results for several m g and m c , shown in Fig. 7, lead to analogous conclusions: scaling described by approximate theories works well, especially for D SS , except for the "resonance" intermediate region of d g . In particular, this scaling can be used to estimate the effect of uncertainty in m g and m c on the final simulated results.
All the results above are for several side scattering signals. In Fig. 8 results for other scattering quantities are shown, namely Q ext and <cosθ >. One can see that these quantities only moderately depend on a particular granule placement, as indicated by relatively small SDs. However, they non-trivially depend on both d g and f. For small d g , Q ext is close to that of a homogeneous sphere with EMT refractive index, which linearly depends on f. Since dependence of Q ext of a sphere on m is oscillating, so is the dependence of Q ext of a granulated sphere on f for small d g . For large d g overall extinction by granules is less dependent on each other, therefore Q ext increases with f for constant d g (if f = 0, Q ext = 1.04). In intermediate region of d g Q ext smoothly changes from one limiting value to another. The asymmetry parameter can be thought of as another measure of side scattering intensity. And one can see We do not discuss other scattering quantities, e.g. angle-resolved Mueller matrix elements, in this manuscript, because it is hard to choose any of universal interest. Any particular application needs its own combination of Mueller matrix elements and range of scattering angles.
As we have seen above, RDG and 2-Born adequately describe the scaling of the side scattering signals. It is much more enlightening to perform a quantitative comparison of these approximate theories with DDA results. We present an example of such comparison in Fig. 9 for the default parameters of the DDA simulations (see Section 2). RDG results are shown for three volume fractions f = 0.02, 0.05, and 0.1, while 2-Born only for the default one ( f = 0.1). One can see that RDG is an accurate approximation for small f, especially for small x g . However, it systematically underestimates I SS for larger f. This discrepancy is due to multiple scattering effects which are significant for larger f and are completely ignored in the framework of RDG. The RDG also has some more pronounced limitations, for example, the RDG result for Q ext do not depend on x g but only on f, contrary to the DDA simulations [see Fig. 8(a)]. Probably, RDG is more accurate for <cosθ >, but it cannot be accurately calculated in the framework described in Section 4 because of the employed assumption of not small θ.
It is important to note that 2-Born results shown in Fig. 9(b) [Eq. (56)] are obtained for zero scattering aperture, while DDA results are for Δθ = Δϕ = 25°. We do not use DDA results for zero aperture because their standard deviations are larger than the values themselves. Since we calculated only 10 different realizations of granule positions for each x g , the errors of mean value are also large. One can see that the agreement between the 2-Born and the DDA is good, especially up to the first maximum. For larger x g the 2-Born systematically underestimates the depolarized intensity, which is due to the neglect of higherorder scattering and the fact that even a single large granule is treated inaccurately in the framework of 2-Born.

Conclusion
We performed extensive DDA simulations of light scattering by a granulated sphere. We calculated total and depolarized side scattering intensity (I SS and I ⊥ ) varying parameters of the model: granule diameter d g , volume fraction f, and refractive index m g , cytoplasm diameter and refractive index and size of the side scattering aperture. Parameters used correspond to the literature data on granulocyte morphology. The general appearance of the maps I SS versus I ⊥ for several d g is the same for different values of other model parameters: two segments with high correlation between I SS and I ⊥ for small and large d g and a connecting region for a narrow range of d g . Both I SS and I ⊥ increase with d g for small d g , have a maximum at d g approximately equal to the wavelength λ, and slowly decay afterwards. Depolarization ratio D SS behaves in a step-wise manner, having almost constant value for small d g and a larger constant value for larger d g .
Analytical expressions, obtained in the framework of the Rayleigh-Debye-Gans and second-order Born approximations, describe these results qualitatively well. Moreover, these expressions quantitatively agree with rigorous DDA simulations for small d g and f. In addition to being extremely fast, these approximate theories give an insight into the light scattering phenomena. In particular, depolarization of granulated sphere is determined by short-range interaction between nearby granules for very small d g , while for d g of order λ and larger it is determined by long-range interactions of all granules. Although, analytical expressions are not accurate enough in many cases, they can be used to construct approximate inversion techniques, which goal is to deduce information about granules from light scattering signals.
We showed that differences between experimentally measured signals of neutrophils and eosinophils [14] can be described solely by the difference in their granule sizes, which agrees with previous phenomenological descriptions of this phenomena. Moreover, calculated D SS quantitatively agrees with experimental results. However, the quantitative comparison is hindered by the uncertainty of m g . Finally, as an example of other scattering quantities we showed and discussed the dependence of extinction efficiency and asymmetry parameter on f and d g . They can also be used to approach the inverse problem.