Mature red blood cells: from optical model to inverse light-scattering problem

: We propose a method for characterization of mature red blood cells (RBCs) morphology, based on measurement of light-scattering patterns (LSPs) of individual RBCs with the scanning flow cytometer and on solution of the inverse light-scattering (ILS) problem for each LSP. We considered a RBC shape model, corresponding to the minimal bending energy of the membrane with isotropic elasticity, and constructed an analytical approximation, which allows rapid simulation of the shape, given the diameter and minimal and maximal thicknesses. The ILS problem was solved by the nearest-neighbor interpolation using a preliminary calculated database of 250,000 theoretical LSPs. For each RBC in blood sample we determined three abovementioned shape characteristics and refractive index, which also allows us to calculate volume, surface area, sphericity index, spontaneous curvature, hemoglobin concentration and content.


Introduction
Determination of such red blood cell (RBC) characteristics as volume, shape, and hemoglobin concentration allows one to assess a blood status and to reveal a number of diseases such as anemia, hereditary spherocytosis, hereditary elliptocytosis, sickle cells disease, and myeloproliferative disorders [1]. Manual observation of cells by light microscopy is widely used for RBC characterization [2], but it neither has sufficient accuracy nor is able to measure hemoglobin concentration. Automatic counting and volume measurement of RBCs is commonly performed with the impedance method [3]. In spite of statistically significant number of measured RBCs, it ignores shape and hemoglobin-concentration differences between measured cells introducing uncertain inaccuracy. A similar technique is the flow cytometry, based on measurement of two-angle light scattering by single cells [4] or angleresolved light scattering [5]. The main limitation of these methods is the need to chemically spherize RBCs prior to their characterization. This requires additional probe preparation, and only the cell volume, but not the shape can be estimated. Both approaches do not verify spherization efficiency that results in uncertainty of determined volume. A number of interference methods do allow measuring volume, shape, and hemoglobin content (refractive index) of individual RBCs [6][7][8]. However, they are based on simplifying assumptions in the light-scattering theory (thus, have hard-to-control accuracy) and are not yet ready for highthroughput applications, such as routine clinical blood analysis.
In this study we modernize the angle-resolved light-scattering approach based on the scanning flow cytometer (SFC) [5,9,10] for complete and rapid characterization of RBCs from light-scattering profiles (LSP). The new method does not require chemical spherization of RBCs and utilizes the solution of the inverse light-scattering (ILS) problem, using the optical model of a mature RBC. We used similar approach for characterization of lymphocytes, platelets, rod-like bacteria, blood microparticles, and milk fat globules [11][12][13][14]. Extension of this approach to the RBCs requires knowledge of their precise shape. There exist a number of simple empiric parametric shape equations, which has been used in optical studies of RBCs [6,15,16]. While most of them has been to some extent verified by fitting to the experimental data, they cannot be used to reliably describe the whole range of mature RBC shapes, since the underlying data is limited to a few hundred cells [6,8]. Therefore, a physically-based model would be preferable. From a mechanical viewpoint, the RBC is a bag filled with hemoglobin, its shape is mostly determined by the minimum of the membrane free energy [17]. It was assumed [18] that the main shape-dependent part of the latter is the membrane bending energy. This idea was further developed [19], introducing the spontaneous curvature to account for the asymmetry between the two parts of the membrane bilayer. This model results in a variety of RBC shapes: cup-shaped, stomatoid, discoid, and dumbbell axisymmetric ones under variation of the spontaneous curvature and the cell volume at constant membrane area. For applications involving many instances of the shape, i.e. for fitting a LSP to solve the ILS problem, the shape model needs to be both realistic and quickly computable from a small number of characteristics.

Optical model of red blood cell
A typical human RBC shape is a biconcave disk with diameter of ~7.5 µm, a thickness of ~2 µm, a volume of ~90 μm 3 , a surface of ~130 μm 2 [6]. We make two assumptions: a) gradient forces near the center of Poiseuille flow in SFC capillary are sufficiently small not to deform a particle, b) contribution of membrane cytoskeleton is important only for highly abnormal shapes [20], which are not considered here. Under these assumptions, RBC tends to acquire the resting shape corresponding to the minimum value of the membrane free energy: with given cell volume V, sphericity index SI = 3(4π) 1/2 VS 3/2 , and membrane spontaneous curvature c 0 , where с p and с m are two principal curvatures of membrane surface, k c is the elastic constant, A is a closed surface. SI is the ratio of the mature RBC volume to the volume of sphered RBC having the same surface area S. The spontaneous curvature describes the curvature of the relaxed membrane (if a small patch of it is cut out) and it can differ from zero due to differences between the inner and outer membrane layers. Note that E(A,c 0 ) = E(kA,c 0 /k), where kA is the homogeneous scaling of surface A by a factor k. Therefore, all dimensionless shapes determined by Eq. (1) can be obtained by varying only SI and c 0 for any fixed volume V, which plays a role of a scaling factor.
Two different numerical methods were used for RBC shape simulation. The first method is an iterative process of energy minimization with simultaneous 3D surface modification by the Surface Evolver software [21]. Optimization process, starting from some initial shape, is realized with constraints on volume, surface, and spontaneous curvature. While Surface Evolver can deal with arbitrary shapes, including cup-shaped RBCs, it requires human supervision, and its iteration process is cumbersome. Therefore, we used Surface Evolver only as a reference to verify the second method, based on the simulation of axisymmetric 2D profile, described by Deuling & Helfrich [19]. Additionally assuming that the RBC shape has mirror plane (xy) perpendicular to the symmetry axis (z), the problem can be transformed into the Euler-Lagrange equation. The latter is solved by the shooting method, in which the boundary conditions are varied to obtain smooth coupling of two branches. Examples of obtained 2D RBC profiles at different values of с 0 and SI are shown in Fig. 1. Here and further x is equivalent to the axial cylindrical coordinate ρ.
Numerical tests show ambiguities in solutions: fixed set of characteristics (V,SI,c 0 ) may correspond to two or more different shapes with minimum energy, see, e.g., Fig. 1(d)). Therefore, we propose another set of characteristics: d, h 1 , h 2diameter, bottleneck, and thickness respectively (Fig. 2). Their one-to-one correspondence to shape was shown by numerical tests with random variation of these characteristics (data not shown). Moreover, h 1 /d and h 2 /d uniquely determine the dimensionless shape, which can be further scaled to any required d. Unfortunately, the new characteristics add additional level of complexity. While finding them for any specified shape is numerically fast, inverse procedure requires finding those (V,SI,c 0 ), for which the shape has required (d,h 1 ,h 2 ), i.e. a two-level optimization.  To alleviate this complexity we use an analytical approximation for the shape profile. To find the latter, which is both convenient and accurate, we considered the following popular parametric equations: Cassini ovals [22], equation of Borovoy et al. [15], model of Fung et al. [6], and proposed an extended Fung model of variable order N: where the free parameters are C 0 -C N1 and d is scaling factor. These shape equations were tested by trying to fit realistic shapes corresponding to the minimum bending energy. The extended Fung model of fifth order provides the best fits over a wide range of (V,SI,c 0 ), even for highly non-spherical particles (V = 90 fl, SI = 0.55, c 0 = 0 μm 1 ) or (V = 90 fl, SI = 0.7, c 0 = 1.2 μm 1 ), while other parametric equations do not fit well the reference profiles (data not shown). Moreover, this model accuracy is within the level of the thermal fluctuations of the RBC shape, which has been measured by common-path optical tomography to be about 40 nm [8]. To further employ the fifth-order extended Fung model we obtain coefficients C 0 -C 4 as functions of h 1 /d and h 2 /d; specifically, we set C 0 = h 1 /d (by its definition), and evaluate C 1 -C 4 by polynomial approximations: where the coefficients , k ij p are determined by a linear regression and tabulated in Data File 1.

Solution of ILS problem
The ILS problem was solved analogously to that described in [14] for characterization of blood platelets. In short, the problem is transformed into the global minimization of the weighted sum of squares: where Q is a vector of model characteristics (Table 1), I th and I exp are theoretical and experimental LSP, respectively, k is the number of LSP points (in the range of θ from 10° to 60° with step of 0.5°), and w(θ) = (1°/θ)exp(2ln 2 (θ/54°)) is the weighting function [11]. Theoretical LSPs were simulated using the discrete dipole approximation, specifically opensource code ADDA v.1.2 [23], using 10 dipoles (volume discretization elements) per the wavelength, as justified previously [24]. Since the calculation of a single I th (θ,Q) takes about 1 minute on a single processor core, the optimization was performed by nearest-neighbor interpolation with a precomputed database of LSPs. The latter contains 250 000 theoretical LSPs with characteristics randomly distributed in the intervals shown in Table 1 [25] at λ = 660 nm. The range of orientation angle Ψ (between the RBC symmetry axis and propagation direction of the incident light) was chosen under the assumption that RBCs are usually oriented with their diameter along the flow direction [9]. The refractive index of the medium (saline) is n 0 = 1.333. Those simulations were performed at the Supercomputing center of the Novosibirsk State University.
Comparing I exp with all I th from the database, we not only find the best-fit characteristics Q 0 that minimizes Φ(Q), but also obtain an approximate description of the whole surface of Φ(Q). The latter is used to calculate probability density function P(Q) over characteristic space for a given experimental LSP through the Bayesian approach. P(Q) is further used to calculate mathematical expectations (ME) μ = <Q> (generally different from Q 0 ) and standard deviations (SD) of characteristic estimates.

Experiment and results
Blood sample were taken from a healthy donor by venipuncture and collected in a vacuum tube containing EDTA as anticoagulant. The sample was 1000-fold diluted in 0.9% saline and measured with the SFC, described in detail elsewhere [9]. The particular SFC was produced by Cytonova LLC (Novosibirsk, Russia) and is equipped with 40 mW laser of 660 nm (LM-660-20-S, Laser Technologies LLC, Russia) for generation of LSP of individual particles. Another 50 mW laser of 488 nm (DL-488-050, CrystaLaser, USA) is used to generate the trigger signal. The measured LSP is expressed as: where S 11 is the Mueller matrix element, and θ and φ are polar and azimuthal scattering angles. For each measured LSP the ILS problem is solved (see Section 3) to estimate d, h 1 , h 2 , n, Ψ (ME) and their uncertainties (SD). The derivative characteristics: V, S, SI, с 0 , hemoglobin concentration and content are also calculated (both ME and SD). Typical experimental LSPs together with best-fit theoretical ones and estimated characteristics are shown in Fig. 3. Note that the quality of the fit and uncertainties of determined characteristics vary greatly between individual RBCs, which can be due to both model errors (the deviation of the real RBC shape from the used model) and instrumental ones, e.g. imperfect trajectory of RBC in the flow [11].

Discussion and conclusions
We developed a method for characterization of mature RBCs, based on measurement of LSPs from single cells in flow with the SFC and subsequent solution of the ILS problem for each LSP. In addition to the best-fit characteristics of each RBC, we estimated errors (uncertainties) of these values. The method features sub-diffraction resolution in RBC sizing coming down to 0.1 μm in some cases. More specifically, the median uncertainties of measured diameter and thickness were 0.25 µm and 0.07 μm, respectively. Such spatial resolution opens a way for precise studies of processes involving RBCs, i.e. lysis, activation, etc. Overall, the method allows measurement of volume, surface area, diameter, maximal and minimal thicknesses, sphericity index, spontaneous curvature, refractive index, hemoglobin concentration, and hemoglobin content of mature RBCs. These 9 characteristics (excluding the refractive index) lead to at least 27 different RBC sample parameters (distributions of sample over the characteristics and corresponding mean values and distribution widths) for hematological survey, whereas modern hematological analyzers provide only 8 RBC parameters. Such detailed characterization of RBCs with high speed and precision can be used in routine clinical diagnostics to improve screening examinations. Moreover, the proposed easy-to-use shape model will be useful in other studies (optical, mechanical, etc.) requiring a realistic description of mature RBCs. Further research can be related to non-discoidal shapes (cup-shaped ones, sickle cells) and to including cytoskeletal proteins in RBC shape model to deal with abnormal cells.