Quantitative discrimination of biological tissues by micro-elastographic measurement using an epi-illumination Mueller matrix microscope

: We propose a method for estimating the stiﬀness of bio-specimens by measuring their linear retardance properties under applied stress. For this purpose, we employ an epi-illumination Mueller matrix microscope and show the procedures for its calibration. We provide experimental results demonstrating how to apply Mueller matrix data to elastography, using chicken liver and chicken heart as biological samples. Finally, we show how the histograms of linear retardance images can be used to distinguish between specimens and quantify the discrimination accuracy.


Introduction
The determination of tissue mechanical properties is of large interest in clinical applications for the diagnosis and detection of many diseases, since these properties vary with the pathological condition of the tissue. When external forces such as tension, compression, or shear are applied, the stress and strain of a pathological tissue is different than healthy tissue. For example, cancerous tumors are often stiffer than normal tissue [1].
Micro-elastography provides information on tissue micro-structure and on tissue behavior under applied stress. It is relevant to diagnosis and identification of tissue disorders and disease-related symptoms and thus can act as a method to distinguish between normal and diseased tissue. Elastography techniques have been used by many research groups [2,3] by employing ultrasonic and MRI methods. Elastography has been combined with optical coherence tomography [4][5][6] to determine mechanical behavior of tissue under a compressive load. The measurement is qualitative and the nature of elastograms maps tissue strain with spatial resolution. The above mentioned methods requires complex hardware, and measurement artifacts can be a problem if the data and image processing are not done with great care. Polarization imaging provides additional features such as structural, biochemical and functional information of the light interacting medium compared to conventional intensity imaging methods and therefore can be helpful for clinical applications. Recently, Buchta et al. [7] proposed a shearing-interferometry-based elastography method for determining the elastic parameters and localized stiffness inhomogeneities of soft tissue. He et al. [8] proposed a method to compare among different tissues by using frequency distribution histograms of the Mueller matrix elements. Du et al. [9] show that by using Mueller matrix polarimetry, characteristic features of cancerous tissues can be differentiated from normal tissue.
Biological tissue is anisotropic where collagen fiber structure is present in tissue, such that fibrous tissue can be differentiated from other tissue types using linear retardance measurement. Since abnormal tissue shows different fiber structural properties than does normal tissue, linear retardance thus can be used as an effective tool for diagnostic purposes [10][11][12][13]. Mueller matrix imaging contains information of retardance, diattenuation, and depolarization allowing for a better understanding of polarization properties of sample and thus a Mueller matrix microscope is capable of quantitative analysis of a tissue's birefringence and its direction in the micro scale. Epi-mode measurement is useful compared to transmission mode for measuring opaque samples such as biological specimens.
The most commonly used instrument for Mueller matrix polarimetry is perhaps the dual rotating retarder technique originally proposed by Azzam [14]. Arteaga et al. [15] introduced a transmission mode Mueller matrix microscope using the dual rotating retarders methodology. For elastography measurement, the specimen should have a thickness of the order of mm and in that case, a transmission Mueller matrix microscope can not be used because of light penetration and therefore it is required to build a Mueller matrix microscope in epi-illumination mode. We introduce a epi-illumination Mueller matrix microscope instrument for the quantitative discrimination of elastographic properties of tissue.
In the discussion below, we survey the measurement model for our epi-illumination Muller matrix microscope system and explain the calibration method for it. In order to demonstrate the effectiveness of our measurements, we record Mueller matrix images of a rubber sample while applying different amounts of stress. After using the rubber sample to verify the accuracy of measuring Young's modulus using the Mueller matrix microscope system, we establish an empirically determined linear relationship between the sample linear retardance in reflection and the applied stress. This relationship is described by what we call the sample's "stress-retardance sensitivity coefficient" and by comparing the stress-retardance sensitivity coefficient of different tissues, we show that it is possible to discriminate tissue types even when they appear same to the eye.

Methodology and calibration of the instrument
Our Mueller matrix microscope is depicted in Fig. 1. The light source used in this experiment is a halogen lamp (Ocean Optics HL-2000-HP) which is transmitted through an optical fiber (Ocean optics, P1000-2-UV/VIS, 1000 µm dia). The light beam is collimated by a collimating lens (Edmund Optics, 43-902, NA= 0.25) and then passed through a bandpass filter (Edmund Optics, 65-098, 12.5 mm dia) having a spectral bandwidth 10 nm, centered at 550 nm. The focal length of the collimating lens is 15.37 mm and the diameter of the collimated beam is 8 mm. The polarization state generator (PSG) is composed of a horizontal polarizer and a rotating quarter wave plate (QWP). The polarization state analyzer (PSA) is similar to the PSG with the components placed in reverse order. The polarizers are Glan Thompson type having 12 mm diameter with an extinction ratio greater than 100000:1 and the waveplates are zeroth-order QWPs. The two retarders in the PSG and PSA are mounted independently with motors to rotate at rotation speeds of θ(t) and 5θ(t) respectively. The beam splitter is a non-polarizing beam splitter used to pass light from the PSG into the microscope objective lens. The microscope objective is a plano infinity-corrected long working distance objective lens (Mitutoyo, 10X, NA=0.28). Imaging is performed using a cooled CCD camera (Bitran BQ-86M, 1360×1024 pixels, 16 bits).
The exposure time of the camera can be varied depending on the requirement of the measurement. Dark current is measured before each measurement and is subtracted from the signals at each pixels during measurement. Averaging of five measurements are performed to improve the signal to noise ratio and to increase the precision of the measurement. Total thirty six images are taken and then these images are used to determine the Mueller matrix elements of the sample. The thirty six measurements are spread over a 180°rotation of the PSG. Without any calibration of the system, the Mueller matrix elements can be found from the Fourier coefficients of the output intensity and is given in detail in Ref. [16].
We demonstrate the calibration of our reflection-mode Mueller matrix polarimeter. Since the retarders used in the proposed configuration are only approximately achromatic in the visible wavelength range, working at different wavelengths requires recalibrating the retardance of the wave plates. Also, the main source of errors in a dual rotating retarder instrument is the azimuthal error in the components. Goldstein et al. [17] proposed a calibration method for a transmission configuration by taking air as a reference sample. In reflection mode, this is not feasible. Instead, we take a reference measurement using a mirror as the sample, and from this calibrate the retardance of the retarders and azimuth angles of the components. We also calibrate the effect of the beam splitter in both reflected and transmission paths through the beam splitter, since it causes an amount of linear retardance and diattenuation. Since the polarization effects due to the microscope objective are very small, therefore they are ignored in the calibration process.

Calibration of retardance error of waveplates and azimuthal error of the components
Let us define the linear retardance magnitude and fast axis orientation angle of the linear retarder as δ and θ respectively. Let us consider retardance errors of the two retarders as, 1 = δ 1 − 90°, 2 = δ 2 − 90°. If we write 3 and 4 are azimuthal errors of the retarders and 5 is azimuthal error of the analyzer, then the output Stokes vector at the camera can be written using Stokes-Mueller calculus as where S in is the input Stokes vector of the incident light, P n (φ) be the Mueller matrix of the nth polarizer having transmission axis oriented at an angle φ, and R n (θ) be the Mueller matrix of the nth linear retarder. The reference sample is taken as mirror at normal incidence whose Mueller matrix is given by Using this theoretical model, the corresponding output intensity is Fourier transformed to determine the Fourier coefficients a n and b n . Therefore, the azimuthal errors of the components can be expressed as and retardance errors can be explicitly written as By using this calibration method, the azimuthal angles and retardance of the waveplates are retrieved and then we adopt the equation given in Ref. [17] to retrieve the Mueller matrix elements.

Calibration of the beam splitter
Although the beam splitter is considered as a non-polarizing element, it has linear retardance and diattenuation which can cause artifacts in the measurement results of the experimental Mueller matrix if not compensated-for. To remove the linear diattenuation and retardance effects due to the beam splitter, calibration of beam splitter is also performed.
In the experimental procedure, the light first reflects from the beam splitter, then interacts and reflects from the sample and finally transmits through the beam splitter. In the whole process, the measured Mueller matrix (M measured ) is a multiplication of the Mueller matrix of the beam splitter for reflection (M BSR ), the Mueller matrix of the sample (M S ), and the Mueller matrix of the beam splitter due to transmission (M BST ), and therefore can be written as Complete Mueller matrices of the beam splitter was measured using a commercial Mueller matrix polarimeter (Axoscan) [18] in both reflection and transmission and are shown in Table  1. The AxoScan system measures polarization properties at a single point (not imaging) of the sample and can measure in the visible to near infrared range (400-800 nm). The measurement error in the Mueller matrix elements by the Axoscan system is 0.1% and the system has a precision of 0.01% for 40 set of measurements at our working wavelength of 550 nm. In transmission configuration, the measured linear retardance is 1.42°and diattenuation is 0.281 whereas in reflection mode the beam splitter has a linear retardance of 2.36°and diattenuation of 0.270.
The measured Mueller matrix of the beam splitter is sensitive to the optical geometry of the commercial polarimeter. However, the same geometry was used during measurement of the beam splitter with the commercial polarimeter as with the Mueller matrix microscope.

Determination of stress-retardance sensitivity coefficient
We next show how the linear retardance values obtained from a sample's Mueller matrix can be used to estimate the sample's stiffness. After passing through a birefringent material, the two components of the electric field vector along the two principal stress directions experience a different refractive index. Using the stress-optic law, the magnitude of linear retardance ∆ expressed in radians can be written as where n x and n y are the refractive indices along the two principal directions, G is the stress optic coefficient expressed in radians/unit pressure, σ x and σ y are the amount of stress applied along the two principal axes expressed in Pa, d is the effective optical path length of light through the sample, and λ is the wavelength of the light source. By applying stress in one principal direction (σ x = σ, σ y = 0), Eq. (7) simplifies to: The strain can be calculated as the change in length (∆L) divided by the original length (L). The applied stress σ and the produced strain are related to Young's modulus E by If we measure the strain of a sample as a function of applied stress, we can fit a line to the data to determine Young's modulus for the material. By substituting Eq. (9) into Eq. (8), we obtain  where we define the "strain-optic modulus" R = 2π G E d/λ. From the slope of the retardancestrain plot fitted from the same experiment but this time from the reflected Mueller matrix data, we can also determine the value R. Since strain is dimensionless, R has the same units as ∆.
We call the coefficient a the "stress-retardance sensitivity coefficient". Since Young's modulus E is expressed in Pa and R is in radians, therefore the unit of a will be Pa/rad.

Instrument validation
In order to experimentally validate the Mueller matrix microscope instrument, we measure the spatially resolved Mueller matrix of a mirror at normal incidence. The Mueller matrix element images and element standard deviation images of the mirror are shown in Fig. 2. All of the Mueller matrix elements are normalized by m 00 . Each standard deviation is obtained from a set of five measurements. The mean measurement errors of the diagonal and nondiagonal elements of the Mueller matrix elements of the mirror are less than 0.8% and 1.1% respectively.  We next measured the Mueller matrix of a rotating Glan-Thompson polarizer. The Mueller matrix is decomposed by the polar decomposition method [19] and the retrieved linear diattenuation magnitude and orientation for the different orientation axes of the polarizer are shown in Fig. 3. The value of linear diattenuation averaged across the face of the sample (circular region in Fig. 3) is measured to be 1.00±0.01 i.e. a perfect polarizer to within the measurement precision of our system. The ±0.01 error in linear diattenuation represents the spatial standard deviation. The orientation of diattenuation follows the azimuthal axis of the linear polarizer and has an error of less than 0.2°.
Next, a waveplate (λ/4 at 632.8 nm) is measured to quantify the instrument's retardance measurement accuracy at wavelength 550 nm. Figure 4 shows the linear retardance magnitude and relative axis orientation measured while rotating the waveplate from 0°to 180°. By taking the spatial average and spatial standard deviation across the face of the sample, the linear retardance of the waveplate is measured as 105°±2.4°, independent of the rotation angle of the waveplate. The estimated orientation of linear retardance follows the waveplate orientation to within 0.5°of error.

Elastographic measurement
The mount depicted in Fig. 5 is used to measure the strain produced in a sample. Stress is applied by attaching different weights to a pulley. One side of the sample is held fixed on the fixed stage which is kept in place while the other side of the sample is attached to a sliding stage. The   direction of the strain is perpendicular to the incident beam.
To validate the accuracy of the stress-strain measurement using this technique, we first performed a measurement on a rubber material having a known Young's modulus of 10.8 MPa. For a given strain, the complete Mueller matrix was simultaneously recorded for the rubber sample and the polar decomposition method used to retrieve the linear retardance. The spatial average linear retardance for strains from 0.000 to 0.048 are shown in Fig. 6. Linear retardance in the rubber arises due to stress birefringence for the applied stress produced in the sample. To quantitatively evaluate linear retardance images (Fig. 6), the spatial histograms of linear retardance and fitted Gaussians for different applied strain are shown in Fig. 7. The mean value of retardance increases linearly with increasing strain (Fig. 8). Figure 8 also shows the stress-strain measurement result. Fitting the measured stress-strain data to a straight line gives an estimated Young's modulus of 10.5±0.12 MPa. The ±0.12 range in Young's modulus represents the standard deviation after five measurements. The value of the stress-optic modulus R for the rubber sample is determined as 220°using Eq. (10). We can estimate the stress-retardance sensitivity coefficient a, after measuring the Young's modulus E and stress-optic modulus R from their ratio (Eq. (11)) as of 47.7. Next, we perform measurement of biological tissues. Samples of raw chicken liver, raw chicken heart, (shown in Fig. 5) and cured ham were each cut in a rectangular shape of 10 mm×20 mm cross-section and 1.2 mm thickness using a long blade knife. We prepared the tissue in rectangular shape rather than free form shape and then applied force in one direction, such that the deformation in the tissue is approximately uniform over the sample and is confirmed by imaging four different parts of the sample. To determine the stress-retardance sensitivity coefficients of chicken heart and liver, we performed stress-strain and Mueller matrix measurements simultaneously by giving an increasing amount of applied stress. When analyzing the Mueller matrix of birefringent turbid media with Mie-sized scatterers acquired in reflection geometry, the Lu-Chipman decomposition may suffer from limitations due to the assumptions required by this method. In this case, one can make use of the extended polar decomposition method instead. [20]. Figure 9 shows the linear retardance images of chicken liver and heart for different strain value. The mean values of linear retardance are plotted with strain in Fig. 10. The non-zero retardance at zero strain is due to the small amount of birefringence presence in the samples due to their fiber structure.
At lower values of applied stress, the measured strain and retardance increase linearly. Between 0.5 and 1.0 of strain, however there is a sudden change in behavior for chicken liver, where the retardance-strain curve exhibits a break from its initial slope to a new lower slope. Above this point, the retardance drops suddenly and then becomes more or less constant with applied strain. This can happen due to the breaking of the internal fiber structure and tissue damage. We fit the stress-strain data and linear retardance-strain data from the first four data points before when the fiber structure breaks. Fitting the measured data to a straight line through the first four data points by least square fitting, we calculate Young's modulus of chicken heart and liver as 0.73 MPa and 0.14 MPa respectively. Thus, chicken heart tissue is five times stiffer than chicken liver tissue.
Using Eq. (10), the value of the stress-optic modulus R for chicken heart and liver are determined as 230°and 51°respectively. To confirm uniform strain over the sample, we also measured Mueller matrix in four different parts of the chicken heart and liver. The spatial histograms of the linear retardance in the four different portions of the sample almost matches with each other and confirm that the stain is uniform over the sample.
After measuring the Young's modulus E and stress-optic modulus R, we can estimate the stress-retardance sensitivity coefficient a from their ratio (Eq. (11)) as of 3.2 and 2.7 kPa/°chicken heart and liver respectively. While Fig. 10 gives the estimated value for the stress-optic modulus R, computing R from material properties alone (i.e. Eq. (10)) requires knowing the "effective optical depth" of the sample in reflection. This is not same as the sample thickness since the light can not penetrate all the way through. In order to estimate d for our samples, we used a microtome to slice the tissues in thicknesses of 5, 10, 15, and 20 µm and compared the retardance measured from these thin slices to retardance measured on bulk samples. Since the 10 µm thin slices gave results comparable to the bulk samples, we use d≈10 µm for the effective optical depth.
We also measure the stress-strain properties and Mueller matrix image of a cured ham sample. The linear retardance image is extracted from Mueller matrix image for different amount of strain. As before, we get retardance images at different values of stress and fit the stress-strain data and retardance-strain data up to the point where the tissue fibers show damage. The slope of the fitted line gives the Young's modulus of ham as 0.59 MPa [21]. Using the same fitting procedure, the value of the stress-optic modulus R is measured to be 170°. From the Young's modulus E and stress-optic modulus R, the stress-retardance sensitivity coefficient a of ham is determined as 3.4 kPa/°.
In order to demonstrate how the stress-retardance sensitivity coefficient can be used to discriminate between two different tissue types, we show the measured linear retardances for chicken heart and liver at the same value of applied strain 0.2 (Fig. 11). From the fitted Gaussian curves, we can see that chicken heart has a higher retardance than chicken liver for a fixed applied stress. A standard choice for discriminating between the two samples is to draw a line where the two curves intersect-in this case at 27.2°. All pixels to the left can be classified as liver tissue, 3   and all to the right as heart tissue, for which we can calculate a classifier specificity of 0.986 and sensitivity of 0.977. A comparison of different physical parameters in terms of Young's modulus, the stress-optic modulus, and the stress-retardance sensitivity coefficient of the samples used in the present study is listed in Table 2. For the biological samples, it is evident from Table 2 that, the properties E and R closely track one another , so that R can be used as a kind of proxy measurement for E. However when the sample has a very different behavior, such as the case for rubber sample versus all of the biological samples, then there is no clear connection between E and R. Then the stress-retardance sensitivity coefficient a can be useful for discrimination between different samples.   Chicken liver 0.14 51 2.7

Conclusion
We have demonstrated a method to determine the elastographic parameters in terms of a stress-optic modulus and a stress-retardance sensitivity coefficient for differentiation of biological samples. An epi-illumination Mueller matrix microscope and its calibration method is demonstrated that is useful for the measurement of opaque samples. Linear retardance images are retrieved from the Mueller matrix elements and are shown for the biological samples including chicken liver and chicken heart. Using the histogram of the linear retardance images, it is shown that we can differentiate between different tissue structure. We quantitatively distinguished between different biological tissues from the stress retardance sensitivity coefficient by correlating Mueller matrix with the mechanical properties of the tissue. While our technique provides for a means to measure thick biological samples, due to its use of reflected light, rather than in thin slices for transmission, accurate quantification requires calibration. Some tissue samples will allow for deeper transmission in this measurement geometry than other samples will, but if their basic transmission and scattering properties are known a priori, then it seems likely that a calibration method can be used to match retardance to stress in each case. We believe that using Mueller matrix polarimetry technique and by applying strain in the tissue, it is possible to distinguish between normal and affected region of the tissue and therefore can be helpful for early detection of many diseases.