Analysis of polarized diffraction images of human red blood cells: a numerical study

: We carried out a systematic study on cross-polarized diffraction image (p-DI) pairs of 3098 mature red blood cells (RBCs) using optical cell models with varied morphology, refractive index (RI), and orientation. The influence of cell rotation on texture features of p-DI pairs characterized by the gray-level co-occurrence matrix (GLCM) algorithm was quantitatively analyzed. Correlations between the transverse diameters of RBCs with different RI values and scattering efficiency ratios of s- and p-polarized light were also investigated. The correlations remain strong even for RBCs with significant orientation variations. In addition, we applied a minimum redundancy maximum relevance (mRMR) algorithm to improve the performance of support vector machine (SVM) classifiers. It was demonstrated that a set of selected GLCM parameters allowed for an efficient solution of classification problems of RBCs based on morphology. For 1598 RBCs with varied shapes corresponding to normal or pathological cases, the accuracy of the SVM based classifications increased from 83.8% to 96.8% with the aid of mRMR. These results indicate the strong potential of p-DI data for rapid and accurate screening examinations of RGC shapes in routine clinical tests. under the terms of the Optica Open Access Publishing Agreement


Introduction
Red blood cells (RBCs) play an essential role in oxygen delivery and carbon dioxide removal in humans.Mature RBCs contain no nuclei and exhibit biconcave shapes with different sizes and a large degree of deformability which is a key factor enabling them to accomplish their physiological functions [1,2].Previous studies have shown that the alteration of human RBC morphology leads to a number of abnormalities, such as iron deficiency, thalassemia, and sickle cell disease [3,4].The shape flexibility of RBCs significantly increases their tendency to aggregate, which affects microcirculation by increasing the viscosity of blood flowing in the capillaries [5].Among the various methods to quantify shape variations and aggregation of RBCs, the measurement of elastically scattered light has received significant research attention [6][7][8][9][10].In this method, the angular distribution of scattered light by single cells is acquired to evaluate cell morphology and discriminate between different cell types, including normal and abnormal RBCs.Compared with complete blood count and manual examination of blood smears, RBCs can be characterized by analyzing the patterns and light distributions of coherent light scattered by RBCs.Quantitative learning of these patterns enables the extraction of the morphological features of the imaged cells.The large signal-to-noise ratio in elastic scattering also makes it possible to rapidly image scattered light signals by flow cytometry or microfluidic methods for the characterization and classification of various cancer cell lines and RBCs accurately and rapidly by their morphological differences or changes [11][12][13][14].To fully characterize the patterns of scattered light by RBCs of different shapes, new modeling approaches that also promote the development of novel methods of measurement for clinical applications have been investigated [15,16].For example, several studies have focused on simulating and measuring the integrated or angle-resolved distribution of scattered light intensity along a 1D curve, yielding limited information for reliably characterizing the morphological features of RBCs.In contrast, imaging the angular distribution of scattered light on a 2D sensor surface can provide much enhanced information to assess the morphology and changes of RBCs for quantitative assays [17,18].
We have developed a rapid method of polarization diffraction imaging flow cytometry (p-DIFC) to image coherent light scattered by single cells [19].Compared with the existing methods of single cell assay, its benefits lie in the speed of cell processing which can detect up to about 3000 cells label-free per minute and the ability to reveal 3D shape differences among imaged RBCs instead of cell counting and sizing only [20,21].In this method, the scattered light is divided into p-and s-polarized components to acquire one pair of cross-polarized diffraction images (p-DI).The coherent nature of the scattered light makes p-DI data sensitive to the topological changes of cells in each of the three dimensions which is very valuable for recognition of diseases related shape changes in RBCs.So, the p-DI pair contains rich information on the intracellular distribution of the refractive index (RI) of the imaged cell, which can be regarded as homomorphic to the cell morphology observed through conventional microscopy methods.Using the p-DI data, one can distinguish cells of different phenotypes by learning and extracting the embedded features of the diffraction pattern related to cell morphology.To quantitatively correlate the 3D intracellular RI distribution to the measured p-DI data, we further developed a set of tools to simulate light scattering by a cell with a realistic optical cell model (OCM) and calculate p-DI pairs by tracking scattered light through the imaging unit of a p-DIFC system [22].The results of our previous investigations, ranging from single spheres to human cancer cell lines as well as human primary T lymphocytes, validated and demonstrated the effectiveness of these tools [12][13][14][22][23]. Compared to cells with intracellular organelles, RBCs present simpler structures in terms of heterogeneity in molecular types and RI distributions.Mature RBCs have no nuclei and are filled with hemoglobin molecules of different variants, so they can be treated as scatterers with a constant RI value properly weighted among different hemoglobin types.Hence, the extension of our tools to obtain calculated p-DIs by OCMs of RBCs can provide important insights into the dependence of diffraction patterns on the shapes and aggregated configurations of RBCs.In this study, we numerically studied the influence of RBC shape, size, and deformability on the image textures of p-DI pairs in the cases of deformed RBCs due to flow pressure and diseases such as thalassemia and anemia.Furthermore, we applied p-DIs to investigate the classification of normal and abnormal RBCs, and the results were satisfactory.The OCMs of RBCs and simulation results may offer a new means for evaluating the properties of RBCs with potential diagnostic capabilities.

Methods
A schematic diagram of our p-DIFC system is shown in Fig. 1(a).The p-DIFC system has two paths for the collection of scattered cross-polarized light.The polarization direction of the 532 nm wavelength incident beam was adjusted with half-wave plates and a polarizing beam splitter and focused at the core fluid of the p-DIFC system with a focused beam diameter of approximately 30 µm.Light scattering occurs when the RBCs flow through the focus of the incident beam because of the mismatched RI between the cell and core fluid.The side scatter from the imaged cell with a scattering angle of 68-112°and an azimuth angle of 158-202°is gathered by a microscope objective aligned along the x-axis and split into s-and p-polarized beams by a polarizing beam splitter.These two beams were imaged by two CCD cameras.By combining discrete dipole approximation (DDA) and ray tracing, we previously developed an accurate method for simulating the light scattering process, which enabled us to obtain calculated p-DI pairs with OCMs to quantify the correlations between cell morphology and the p-DI textures.Briefly, we first built an OCM for a given 3D RBC structure and RI value for hemoglobin molecules inside RBCs.Then, we used the open-source ADDA simulation software to solve the spatial distribution of scattered light based on the DDA model [24].The output data from ADDA are saved in the form of Mueller matrices for selected scattering directions, and different matrix elements can be linearly combined based on the polarization states of the incident and scattered light.The calculated p-DI pairs were finally obtained using the ray-tracing software Zemax (2007, Zemax LLC) to project the scattered light through the p-DIFC imaging unit.For different RBCs, we only need to update the respective OCMs with a preset cell orientation.shows the cross sections of three types of RBCs with different shapes and the size parameters, including the transverse diameter as d, center thickness as b, maximum thickness as h, and radius for the maximum thickness as c [25,26].The surface of a biconcave RBC can be described by the following equation [27]: where r is the RBC radius in cylindrical coordinates given by √︁ x 2 + y 2 and the coefficients S, P, Q, and R are determined from the values of the size parameters shown in Fig. 1(b).According to previous studies, thalassemic RBCs have smaller values of cell thickness, transverse diameter, and RI than normal RBCs.In such cases, the center thickness b can serve as the most useful feature for quantifying the abnormal morphological characteristics of thalassemic RBCs, which tend to have very small values of b approaching zero, as shown in Fig. 1(b).To describe the thalassemic RBCs using Eq.(1), we define the following dimensionless ratios as the relative morphology parameters for this study: h d = h/d as the relative h thickness, b d = b/d as the relative center thickness, and c d = 2c/d as the relative h diameter.Using the above definitions, the coefficients defined in Eq. ( 1) can be expressed as Note that all the coefficients of Eq. ( 1) are well defined, even for b=0.The parameters h d , b d c d , and d were chosen to best fit the specific cells that were imaged in [4] using tomographic diffractive microscopy.In order to develop an optical model of normally deformed RBCs caused by the pressure drop across the RBCs and the resistance of the capillary wall, we used the deformed axisymmetric model proposed by Zarda for different values of a dimensionless pressure drop ∆P [28].The cross-section is generated by curve warping on the basis of the biconcave RBC under the restrictions of volume invariance, as shown in Fig. 1(b).By rotating it around the axis of symmetry, we can build the 3D structure of a normally deformed RBC.
Figure 1(c) shows the calculated p-DIs of three different RBC models with the polarizations of the incident and scattered light marked.In view of the dark current noise of CCD is less than 0.2% of full scale in sensor's 12-bit pixel depth, the effect of CCD noise was ignored in simulations and subsequent analysis.In this study, we built 3098 RBC models according to the variation in the parameters and orientation of the RBC.According to our previous numerical studies, texture features of p-DIs can vary significantly if the range of angles for rotation of imaged cells is larger than 10 degrees or more [29].But we have observed from the measured p-DI data that cell types of high morphological similarity can be accurately classified.Thus, one must expect the imaged cells in the core fluid of our flow system must be aligned sufficiently in their orientations along a preferred direction [14].Therefore, we defined the orientation as θ by the rotation axis of symmetry of the RBC and the flow direction and limited it to a small range.If the rotation axis of symmetry of the RBC rotates toward the flow direction in an anticlockwise direction, θ is positive; otherwise, θ is negative, as shown in Fig. 1(a).The parameters h d and c d were set to constant values of 0.295 and 0.65, respectively, and d was varied from 4µm to 11µm [30,31].According to [25], the norm-RBC has a diameter of approximately 6-8 µm while macro-RBC has a diameter greater than 8µm and micro-RBC has a diameter smaller than 6µm.Among these models, 2130 RBC models have the same geometrical shape as biconcave RBCs.The parameter b of the thalassemic RBC models varies from 0.18 to 0 with a step size of 0.02.The deformed RBC models were constructed on the basis of biconcave RBCs according to four different values of ∆P: 1.2, 3, 10, and 26.All RBC models were treated as homogeneous, with a constant RI within the cells.Considering the difference in RI of deoxygenated and oxygenated hemoglobin, all RBC models have two RI values of 1.4 and 1.385 for visible light, except for the thalassemic RBC models, and the RBCs are assumed to suspend in blood plasma with an RI value of 1.385 [32,33].Table 1 summarizes the morphology and RI parameters and the total number of RBC models used in this study.
To quantify the textures of the p-DIs of these groups of RBC models, a gray-level co-occurrence matrix (GLCM) was adopted.Its output matrix elements express the relative frequencies of the paired-pixel values along one of the four array directions.There are 15 parameters whose definitions and symbols are given in the Supplement 1 that can be extracted to characterize the textures of each p-DI.Additionally, we calculated the mean value of the pixel intensity (MEA) of the s-and p-polarization diffraction images to characterize the scattering efficiency in different polarization states.Therefore, every RBC has 32 diffraction parameters for quantification analysis.

Analysis of RBCs with biconcave shapes
With the OCMs of biconcave RBCs constructed with the parameters listed in Table 1, we modeled the light scattering by single RBCs using ADDA with 532nm wavelength incident light.The angle-resolved distributions of the scattered light were produced by ADDA in terms of the Mueller matrix.The elements of the Mueller matrix were combined linearly as the light scattered by an RBC to take into account the linear polarizations of the incident and scattered light collected by the imaging unit of the p-DIFC system.By projecting the s-and p-polarized scattered light into the imaging unit and ray tracing, the calculated p-DI pair can be obtained.We further normalized each of the two images in a p-DI pair and combined them into a 2-channel false-color image with s-polarized DI in red and p-polarized images in the green channel.Figure 2 presents selected examples of the combined images for single RBCs of different sizes and RI values with variations in its orientation.
From these images, we can see that the diffraction patterns embedded in a p-DI pair presented as a combined image vary with the parameters of the RBC models and the orientation of the cell model relative to the direction of the incident beam.To analyze the variations of diffraction patterns quantitatively and investigate the correlations between the parameters of the RBCs and their diffraction patterns, we first used SVM with the cubic kernel function to classify these RBC models based on RI in feature space with GLCM parameters, and the results are shown in Fig. 3.The true positive rate shown in Fig. 3 is defined as the ratio between the number of correctly identified RBCs and the total number of RBCs.We can see that the true positive rate is as high as 98.6%.This demonstrates that the p-DIs are incredibly sensitive to the change in RI and the designed classifier for recognizing RBCs with different RI, SVM combined with GLCM parameters, has high recognition performance.
We then evaluated the influence of rotation on the diffraction characteristics of RBCs.From the Fig. 2(c) and Fig. 2(d) we found that the variation of θ affects not only the texture characteristics of p-DIs but also the scattering efficiency of scattered light with different polarization directions.A parameter of fluctuation degree of ratio (FDR) was defined to quantify the variation of s-and p-polarization scattering efficiency ratio given by s-MEA/p-MEA with θ as the angle of RBC orientation relative to the default value of θ=0°as follows Figure 4 shows that FDR varies from -8.4% to 6.1% for θ in a range of -7°to 7°with a positive slope of change.When the θ lies between -5°and 6°the maximum absolute value of FDR is about 5.4%.In addition, Fig. 4 also shows that the effect of θ variation on FDR is similar for RBCs of different transverse diameters.With the nonlinear least squares fit between the scattering efficiency ratio and transverse diameter d of RBCs, the power law relationship between them at different RI values was obtained and shown in Fig. 5 for three ranges of θ.For the smaller range of -5°≤ θ ≤ 5°, the results exhibit a smooth power dependence on d that can be fitted by Eq. ( 4) and Eq. ( 5) below for RBCs with RI values of 1.385 and 1.4, respectively, at θ=0°r atio 1.385 = 6.524d −0.3602 + 0.4252 (4) Fig. 4. the variation of s-and p-polarization scattering efficiency ratio with θ relative to that of θ=0°(all RBC bi with the same RI value of 1.4).ratio 1.4 = 6.635d −0.3365 − 0.7904 The R-square values of both fitting exceed 0.999 and the 95% confidence boundaries of the fitting and the variation in the scattering efficiency ratio caused by θ are also shown in Fig. 5. Once the orientation angle of RBCs increases to the range from -7°to 7°, the power fitting deteriorates but does not affect classification accuracy of RBCs with different RI values.In addition, Fig. 5 shows that the variation in orientation of the RBCs affects the fitting parameters, particularly when the absolute values of θ are greater than 3°.When the absolute value of θ is less than 3°, the scattering efficiency ratio variation is within the 95% prediction interval of the fitting model.Therefore, when the RBC orientation variation is small, the fitting precision and reliability of the fitting are satisfactory.As mentioned above, the texture feature research of p-DI can serve as an efficient method for RBC size and RI distribution analyses.

Classification of RBCs with normal and abnormal shapes
Usually, the examination of abnormal RBC shapes is required to confirm the diagnosis of an abnormality.Therefore, the influence of RBC shape changes on the features of p-DI was studied, as shown in Fig. 6.The selected examples consist of normal shapes, including biconcave and deformed RBCs, and abnormal shapes, including thalassemic RBCs.It shows that the variation in ∆P has a significant influence on the diffraction patterns.The diffraction patterns of RBCs with different b show a high level of similarity to each other.First, we used the same method mentioned above to classify these RBCs based on shapes in the GLCM parameter space with either 16 selected parameters or all of the 32 parameters.From Fig. 7 one can see that the average accuracy of classification with the set of 32 parameters increases by 10.4% and 6.2% respectively in comparison to that with the set of 16 parameters of only p-or s-polarization.But the accuracy values were still relatively low with the averaged value at 83.8% and a standard deviation value at 12.5%.The main false-negative rate was associated with the classification of RBC Sbi and RBC St .Further analysis showed that 189 out of 201 misidentified RBC Sbi models had an RI value of 1.385.The most likely reason is that the p-DIs are the result of the variations in the multiple parameters of the RBC models when we tried to classify the RBCs based on shapes in the GLCM parameter space with all 32 parameters.We introduced the effects of RI, which played a negative role in improving the accuracy rate of target classification.In addition, the continuous and uniform change of b d and the high sensitivity of GLCM parameters to the change in RI increased this negative impact.Therefore, we removed some of the GLCM parameters that were important for RI recognition from the SVM classifier with the aid of the minimum redundancy maximum relevance (mRMR) algorithm [34].The mRMR algorithm quantifies the redundancy and relevance using the mutual information of variables, pairwise mutual information of features, and mutual information of a feature and the response.It ranks all features for their importance and returns scores.A large score indicates that the corresponding predictor is important.This can help us to determine which GLCM parameters should be removed to improve the performance of the SVM classifier based on the shapes.Figure 8 shows the predictor rank and predictor importance scores using RI as the response.This shows that the drop in score between the fifth and sixth most important predictors is large, while the drops after the sixth predictor are relatively small.The small drops indicate that the difference in predictor importance is not significant.Therefore, we removed the top five scores listed in Fig. 8 from the SVM classifier to improve its shape-recognition performance.Figure 9 shows the confusion matrix of the SVM classifiers of RBC Sbi , RBC Nd , and RBC St with the remaining 27 GLCM parameters.We can see that the average accuracy rate increased to 96.8%, and the standard deviation of the accuracy rate dropped to 1.3%.This demonstrated that with the help of mRMR, we can find an optimal set of features of p-DI pairs to identify a specific characteristic of RBCs.

Conclusions
By quantitative and qualitative analyses, p-DI pair-based assay of RBCs allows accurate recognition of 3D shape changes.This study verified that the characteristics of p-DI pairs are related to many RBC morphology and optical parameters, in which the change in RI was an extremely important influencing factor and the rotation of RBCs has similar effect on the diffraction characteristic of RBCs of different transverse diameters.But it is possible that the effect of latter can be significantly reduced for RBCs flowing in a preferred direction by pressure difference.Under different RIs, precise and reliable fitting of the transverse diameter of RBCs based on scattering efficiency ratios of s-to p-polarization were also achieved.Furthermore, by utilizing the mRMR algorithm, 27 GLCM parameters were selected to obtain highly accurate classification results for RBCs against different shapes.Overall, the results of this work can help develop label-free and automated methods for accurate assays of 3D shapes of RBCs.

Fig. 1 .
Fig. 1.(a) Polarization diffraction imaging flow cytometry system schematic (L: laser; WP: half wave plate; PBS: polarizing beam splitter; FL: focusing lens; FC: flow chamber; OBJ: objective).(b) Cross-section and characteristics of the RBC models (bi: biconcave RBC; t: thalassemic RBC; d: deformed RBC).(c) Three calculated p-DI pairs by RBC models with the polarizations of incident beam and scattered light marked.

Figure 1 (
Figure 1(b)shows the cross sections of three types of RBCs with different shapes and the size parameters, including the transverse diameter as d, center thickness as b, maximum thickness as h, and radius for the maximum thickness as c[25,26].The surface of a biconcave RBC can be described by the following equation[27]:

Fig. 2 .
Fig. 2. Examples of normalized p-DI calculated with biconcave RBC models.Each p-DI is labelled with the structure variable on the top left corner.(a) All with the same RI value of 1.4 and a single variable d; (b) All with the same RI value of 1.385 and a single variable d.(c) and (d) All with the same RI value of 1.4, same d value of 7.5µm, and a single variable θ.

Fig. 5 .
Fig. 5. Power fit of ratio (s-MEA/p-MEA) and d (transverse diameter) with dotted lines and filled areas indicating the 95% prediction interval.

Fig. 6 .
Fig. 6.Examples of normalized p-DI calculated with selected RBC models.Each p-DI is labelled with the structure variable on the top left corner.(a) RBC Sbi all with the same RI value of 1.385 and a single variable d.(b) RBC St all with the same d value of 5.5µm and a single variable b d .(c) RBC St all with the same b d value of 0.1 and a single variable d.(d) RBC Nd all with the same RI value of 1.4, the same d value of 6µm and a single variable ∆P.

Fig. 7 .
Fig. 7. Confusion matrix of SVM classifiers of RBC Sbi , RBC Nd and RBC St .(a) with all 16 GLCM parameters of p-polarization.(b) with all 16 GLCM parameters of s-polarization.(c) with all 32 GLCM parameters of p-and s-polarization.

Fig. 8 .
Fig. 8.The rank of predictors by importance score using RI as response.

Table 1 . Morphology and RI parameters and the total number of the RBC models.
Subscripts of groups ID represent the range of diameters and the geometrical shape of RBC: S = micro-RBC, N = norm-RBC, L = macro-RBC, bi = biconcave RBC, t = thalassemic RBC, d = deformed RBC.
a b RI values of RBC at λ=532 nm.c Diameter of RBC without pressure drops.