Determination of collagen fiber orientation in histological slides using Mueller microscopy and validation by second harmonic generation imaging.

: We studied the azimuthal orientations of collagen ﬁbers in histological slides of uterine cervical tissue by two different microscopy techniques, namely Mueller polarimetry (MP) and Second Harmonic Generation (SHG). SHG provides direct visualization of the ﬁbers with high speciﬁcity, which orientations is then obtained by suitable image processing. MP provides images of retardation (among other polarimetric parameters) due to the optical anisotropy of the ﬁbers, which is enhanced by Picrosirius Red staining. The ﬁber orientations are then assumed to be those of the retardation slow axes. The two methods, though fully different from each other, provide quite similar maps of average ﬁber orientations. Overall, our results conﬁrm that MP


Introduction
Collagen is the major component of extracellular matrix (ECM), which determines the mechanical properties of tissues, such as their stiffness and flexibility, among many others [1]. Alterations of the collagen spatial structure are responsible of a variety of serious pathologies, such as scurvy, osteogenesis imperfecta (also called "brittle bone disease"), characterized by extreme bone fragility, or the Ehler-Danlos syndrome involving hyperelasticity of tissues [1]. Accurate characterization of these structures may thus be of great interest to understand the properties of healthy and diseased tissues and possibly to provide innovative diagnostic tools. Very extensive literature is available about collagen structure [1][2][3]. The elementary building blocks of this ubiquitous protein are polypteptidic assemblies also called α-chains consisting of series of triplets [Glycine-X-Y], where X and Y are mostly proline or hydroxyproline, or may be any aminoacids. These chains are synthetized by specialized cells, where they form left-handed helical structures also called triple helices comprising three interlaced α-chains linked by hydrogen bonds. The diameter and pitch of these helices are 1.5 nm and 8.6 nm respectively, for 300 nm typical lengths (type I collagen). After their intracellular synthesis, the triple helices are excreted into the ECM, where their terminal peptides are removed, triggering their polymerization, which eventually leads to the formation of a variety of supramolecular structures. Nowadays, twenty-eight different types of this protein have been identified. Though, in the following, we will consider only type I collagen, which exhibits a fibrillar structure and is by far the most abundant in vertebrates (it accounts for 90 % of total collagen weight). The triple helices get assembled into 10 to 300 nm diameter fibrils, which usually form collagen fibers with typical diameters between 0.5 and 5 µm. The spatial organization of these fibers can take very different forms, to comply with the requirements imposed by the physiological functions of each tissue. For example, the fibers are uniformly oriented in tendons, they form helicoidal structures in bones, while they are highly entangled in apparently random patterns in derm and viscera. Alternatively, collagen I fibrils form superimposed compact "slabs" in cornea, to simultaneously achieve high optical transparency and mechanical stiffness.
Many imaging techniques are available to characterize collagen at various spatial resolutions. At the nanometer scale, the gold standard technique is transmission electron microscopy (TEM) [4,5]. As this technique requires demanding sample preparation procedures, atomic force microscopy (AFM) [6, 7] has been employed as an alternative technique that can be operated in liquid media and/or on living samples [8]. This technique is also suitable for the study of the mechanical properties of collagen fibrils [9]. At larger scales, optical microscopy is particularly well suited, due to its resolution ranging from a few hundred nanometers to a few micrometers, and typical fields of view of several hundred micrometers. Most often, transmittedlight microscopy is used to visualize histological slides stained with Picrosirius red, Masson's Trichrome or other dyes highlighting collagen. In addition to this classical technique, a variety of optical techniques are now available, such as reflectance confocal microscopy [10] and optical coherent tomography (OCT) [11], which can provide 3D images of unstained ex vivo and in vivo samples. Alternatively, Second Harmonic Generation (SHG) microscopy [12] is a multiphoton technique known to provide 3D visualization of unstained collagen fibers with unequalled specificity and contrast [13][14][15][16][17]. This method is currently the "gold standard" for in situ studies of unstained collagen tissues. Otherwise, polarized microscopy, coupled with suitable tissue staining by Picrosirius red, has been shown to be quite effective for the visualization of fibrillar collagen in histological slides [18,19]. Polarized microscopy is almost always implemented by simply inserting in a standard microscope two linear crossed polarizers, or two circular polarizers with opposite handnesses, one in the illumination beam, and the other after the sample on the detection path. Collagen fibers are then easily detected as they appear shiny while collagen-free regions remain dark. Nevertheless, this "simple" polarized microscopy cannot easily provide the distribution of orientations of collagen fibers in the plane of the slide. On the one hand, circular polarizers reveal the presence of linear retardation but do not provide any information about the azimuth of the slow axis; on the other hand, with orthogonal linear polarizers, the observed contrasts depend on both the retardation and the azimuth of the slow axis. These parameters can be disentangled only by combining images taken at various azimuths of the polarizers. Moreover, these rules apply only if retardation is the only polarimetric effect to be considered, which is not the case with Picrosirius Red staining, which introduces significant diattenuation.
In contrast to these simpler techniques, Mueller polarimetry (MP) can fully characterize the polarimetric response of any sample [20,21]. The potential of this technique for biomedical diagnostics is currently being actively investigated in various conditions (in vivo, with fresh or fixed tissues, phantoms, histological slides...) [22][23][24][25]. In principle, the azimuth of the slow axis (as well as all other polarimetric parameters) can be readily obtained by using MP imaging. The purpose of this work is to validate MP for the orientational characterization of collagen in histological slides as compared to SHG, the current "gold standard" for collagen visualization.
This paper is organized as follows. In section 2, we present the theoretical background of Mueller polarimetry. In section 3, we briefly describe the setups used for SHG and Mueller microscopies, and we highlight the image processing procedures developped to extract the fiber orientations from the MP and SHG images. The results are presented and discussed in Section 4. Section 5 concludes the paper.

Theoretical background
The Mueller matrix (M) is the transmission matrix of the Stokes vector (S) that fully describes the light polarization state: where the Stokes vector is given by: with I the total light intensity and I x , I y , I +45 , I −45 , I LC , I RC the intensities measured through an analyzer oriented along the x, y, +45 • and −45 • directions or transmitting left (LC) or right (RC) handed circular polarization. The Stokes-Mueller formalism is needed to properly describe partially polarized states and/or depolarizing systems. The degree of polarization ρ of a partially polarized state described by the Stokes vector S defined in Eq. (2) is given by: and varies from 0, for totally depolarized states, to 1, for totally polarized states. To keep this paper reasonably self contained, we first recall the essential features of elementary polarimetric properties of any sample. • For diattenuators the transmission varies with the incident polarization between two extremal values, T min and T max . These extrema are reached either for two orthogonal linear polarizations or two counter-rotating circular polarizations. The scalar diattenuation D is defined as: implying that 0 < D < 1.
As circular diattenuation is usually extremely small in biological tissues (at least in the visible) [23], in the following we will consider only linear diattenuators, and will call θ D the azimuth of the linear polarization corresponding to T min , which is also called the low-transmission axis. Then the diattenuator Mueller matrix is: where s = sin 2θ D , c = cos 2θ D and m 11 = 1 2 (T max + T min ) is the transmission for incident unpolarized light. This form is given, among many other authors, by Lu and Chimpan [26] who, however, take the high transmission axis as the reference axis, resulting in opposite signs for s and c.
• Retarders introduce a phase shift R between two orthogonal linear polarizations, or two counter-rotating circular polarizations. This phase shift is expressed in radians or degrees. Again, we will consider only linear retardation, which is largely dominant in biological tissues [23], and we will call θ R the azimuth of the polarization giving the maximum phase retardation, i.e. the slow axis. The explicit form of the corresponding Mueller matrix, is again given in [26]: where s = sin 2θ R , c = cos 2θ R and m 11 , the transmission for unpolarized light, is independent of the retardation parameters. As Lu and Chipman take as the reference axis the fast axis, in Eq. (6) the signs of s and c have been changed with respect to [26].
• Depolarizers reduce the degree of polarization of the incident light. In other words, they introduce some "disorder" in the motion of the electric field in the plane perpendicular to the direction of propagation. Several processes can lead to depolarization, among which multiple scattering is dominant in biological tissues. In the following, we will consider only this type of depolarizers, with the additional assumption that the scattering is isotropic. For such depolarizers, the Mueller matrix takes the very simple form: where a and b respectively describe the depolarization for linearly and circularly polarized incident light. The absolute values of both a and b are less or equal to 1. The smaller these absolute values, the larger the depolarization.
In the most general case, it is not possible to attribute unique properties of diattenuation, retardance and depolarization to a given Mueller matrix. Effective parameters describing these properties can be retrieved by various procedures, among which the most popular is the serial matrix decomposition proposed by Lu and Chipman [26] : where M ∆ , M R , and M D are respectively the Mueller matrices of a pure depolarizer, a retarder and a diattenuator. The polarimetric parameters of these elementary components are then attributed to the input matrix M. However, six other decompositions analogous to that described by Eq. 8 can be implemented as well, by changing the order of M ∆ , M R , and M D in the matrix product [27]. Each of these decompositions provides its own values of the depolarization, retardation and diattenuation parameters. Generally, this ambiguity cannot be overcame, except when the system actually consists of a diattenuator, a retarder and a depolarizer traversed in a well defined order. Of course, in practice this is seldom the case, and a "fundamental" ambiguity remains in the determination of the polarimetric parameters to be attributed to the most general input Mueller matrix M.
Fortunately, for thin samples such as histological slides, depolarization is quite small, due to the absence of multiple scattering. Moreover, the Mueller matrices of such samples are quite close to identity matrices, as diattenuation, retardation and depolarization are also quite small. In this situation, D << 1 and R << 1 rad. Thus in Eqs. (5) and (6), we can neglect the quadratic terms in D and R and keep only the linear ones. Equations (5) and (6) then reduce to: and while the Mueller matrix of isotropic depolarizers is obtained by setting a = 1−ε l and b = 1−ε c with ε l , ε c << 1. Finally, the Mueller matrices of weak diattenuators, retarders and depolarizers commute and their product is given by: where each matrix element corresponds to a specific physical property [28]. The ambiguity mentioned above for the general case is thus removed for Mueller matrices close to identity. In practice, if an experimental Mueller matrix is close to identity, all polarimetric parameters are readily obtained from its elements, provided this matrix actually exhibits the transposition properties appearing in Eq. (11), namely invariance for the first line and column, and sign reversal for the other non diagonal terms.

Sample preparation
Uterine cervical cones were collected from patients and processed to obtain histological slides using usual standard protocol. Briefly, tissues were fixed in neutral formalin for 24-48 hours and embedded in paraffin. 10 µm-thick sections were cut, deposited on a glass slide and stained with Picrosirius red during 5 min. Note that we used a glass coverslip to avoid birefringence induced by usual plastic films. Two different slides from two different patients were visualized sequentially by use of MP and SHG microscopies. Black ink spots on the histological slide combined with specific patterns in the images were used to visualize the very same region of interests (ROI) by the two methods. We verified that Picrosirius Red staining did not affect SHG imaging that is performed usually in unstained tissues.

Instrumentation
The principles of SHG and MP microscopies are illustrated in Fig. 1. In the SHG setup (Fig  1(a)), the primary source is a mode-locked Titanium Sapphire laser operating at 860 nm, and delivering 150 fs duration pulses at 76 MHz repetition rate. The actual excitation power is controlled by means of a halfwave plate mounted on a motorized rotation stage and a Glan Laser polarizer. We use typically 15 mW at the objective focus in our experiments. The beam then passes through a XY scanning system comprising two mirrors set on galvanometric mounts. The mirrors are imaged onto the rear pupil of a 20X, 0.95 NA water immersion objective in order to illuminate the objective entrance pupil over its full aperture. Due to its quadratic dependence on the local intensity, the SHG signal is created essentially at the focal point. This focal point scans the sample laterally when the incident beam is scanned angularly, thus providing a 2D image of the sample in the objective focal plane. The objective can also be moved in the axial direction to provide stacks of tomographic images allowing full reconstruction of the fiber 3D orientations.
When illuminated with IR short pulses, some specific tissue components generate two photon excited fluorescence (2PEF), which is the most common source of contrast in multiphoton microscopy [29]. SHG and 2PEF are easily separated by suitable spectral filters, as 2PEF occurs at larger wavelengths (typically around 500-600 nm vs 430 nm for the SHG). In practice, SHG is highly specific of fibrillar collagen, which lacks any center of symmetry, a mandatory condition for second harmonic generation. In this work, we mostly use the SHG signals for direct visualization of the collagen fibers, while the 2PEF signals are used for visualization of the tissue physiology when necessary. To make sure that the SHG signal is independent of the fibers azimuthal orientations (in the xy plane), the IR beam polarization is converted from linear to circular by a quarter wave plate prior to entering the objective. For practical convenience, in our experiment, 2PEF signals are collected in the backward direction via dichroic mirrors, while SHG signals are collected in the forward direction via a condensor. Both signals are then detected by photomultiplier electron tubes, in a photon counting mode, with 10 µs integration time on each pixel, resulting in 9 s typical acquisition time for our 960 x 960 image.
This setup is described in more details in [30]. The MP microscopy setup is outlined in Fig. 1(b). This setup is basically a standard transmission microscope, illuminated by a classical white source (halogen lamp) to allow full field imaging by means of a CCD camera at the output. The wavelength of operation is selected by an interference filter (typically 40 nm FWHM) placed just before the CCD.
The illumination beam polarization is controlled by a Polarization State Generator (PSG), which sequentially generates four polarization states described by four linearly independent The PSG comprises a linear polarizer, two ferroelectric liquid crystal modulators (FLCs) and a true zero order quarter wave plate at 633 nm set between the two FLCs. The PSA is a "mirror image" of the PSG, with the same elements in reverse order. The FLCs are optically equivalent to retardation plates whose orientations can be switched by an electrical signal between two values, separated by 45 • . Then the PSG (and the PSA as well) generates its four polarization basis states by sequentially switching the FLCs over the four possible sets of orientations. The quarter wave plate is added to make these states "as linearly independent as possible" to optimize the noise propagation from the raw intensity measurements to the final Mueller matrices. More details about the Mueller microscope and on the PSG and PSA based on ferroelectric liquid cystals can be found in [25, 31] . Figure 2 shows a typical example of Mueller image of an histological slide of cervical tissue. All elements m i j but m 11 are normalized by m 11 so that the polarimetric response is evaluated independently of the overall transmission of the slide, which is given by the m 11 matrix element. In the following, the normalized elements will be called m * i j . The diagonal elements are offscale, but the information about the (very weak) depolarization carried by these elements is not relevant for the determination of collagen fiber orientation. We thus focus on the off diagonal elements. The observed trends are the following :

Mueller matrix image processing
1. m * 12 = m * 21 take their largest positive (resp. negative) values on "filaments" oriented horizontally (resp. vertically). As a first result, the simple transposition properties predicted by Eq. (11) are obeyed (see Fig. 2), so that we can extract the polarimetric parameters R, D, θ D and θ R without ambiguity. The next issue is then to relate these optical parameters to the collagen organization. Given that the Sirius Red molecules are linear, with conjugate bonds along their molecular axis, these molecules get bonded and aligned with collagen fibers, if any. Accordingly, Picrosirius Red has been reported to enhance the birefringence of collagen fibers, which shows that collagen fibers and Sirius Red molecules have the same retardance axes [18,19]. Considering the physics of retardance process, we then expect that the retardance slow axis (high refraction index) is along the conjugate bonds, that is along the Sirius Red molecular axis. The same rationale applies to the diattenuation: Sirius Red molecules exhibit a stronger absorption along the conjugate bonds, which means that the low transmission axis is collinear to the molecular axis. Consequently, the stained collagen fibers are expected to produce a diattenuation D and a retardance R whose low transmission and slow axes are collinear, and aligned along the collagen fibers. It means that, for a fiber whose projection is oriented at an angle θ (counted counterclockwise from horizontal), the azimuths θ D and θ R of the corresponding diattenuation and retardance are expected to be both equal to θ . Then the angular dependences on cos 2θ defined in Eq. (11) for the elements m * 12 , m * 21 , m * 34 and m * 43 implies that the extremal values will be reached for horizontal and vertical fibers, while for m * 13 , m * 31 , m * 24 and m * 42 the dependence in sin 2θ locates the extrema at ±45 • . The trends observed in Fig. 2 clearly support this hypothesis. The similarity of the patterns appearing in m * 12 and m * 43 on the one hand, and in m * 13 and m * 24 on the other hand, together with the vanishingly small values of m * 23 and m * 32 provides clear evidence that the diattenuation and birefringence are linear, with coinciding low transmission and slow axes. Thus the "filaments" appearing in Fig. 2 can safely be attributed to fibers (or fiber bundles), and their azimuthal orientation in the plane of the slide can readily be determined from the Mueller image. However, to get reliable values of the fiber azimuthal orientation, the scalar diattenuation D and/or retardance R must exceed the noise level measured on a "blank" slide (slide without any tissue). A typical result regarding retardance is shown in Fig. 3. The left panel Fig. 3(a) shows an histogram of the m * 24 element taken on a blank sample. As expected, this histogram is gaussian, with zero mean value. The middle panel Fig. 3(b) shows the scalar retardance R evaluated from the same blank sample. The corresponding histogram is no longer a zero mean one, with an upper limit of approximately 2 • . This shape is easily understood if one keeps in mind that:

m *
which implies that a zero-mean noise actually "pushes" R towards positive values. Finally, the right panel (Fig. 3(c)) shows the histogram of R for the example shown in Fig.  2 with values up to 10 • . Taking into account the results of the blank sample, we consider that pixels for which R < 2 • correspond to background noise and we evaluate the fiber orientation only for the pixels for which R > 2 • . Figure 4 shows the images of the scalar retardation and of the azimuthal orientation of the slow axis for the whole Mueller image, using this filter (R > 2 • ). As the slow axis can be taken in both directions, we used a RGBR scale from 0 • to 180 • (again, counted counterclockwise from the horizontal). Figure 4(b) is then expected to provide the azimuthal direction of collagen fibers within the histological slide.

SHG image processing
Several image processing methods have been developped to extract the orientation of collagen fibers from SHG images [32][33][34][35]. In this work, we implemented a specific procedure as illustrated in Fig. 5. The raw SHG image of the slide (Fig. 5(a)) is obtained as a maximum intensity projection of the SHG image stack acquired every 1µm within the slide, with 0.5 µm lateral pixel size. First, SHG images are filtered using a mean filter (radius 2 pixels) to remove most of the background noise. Second, a morphological opening [36] is applied using a structuring element (Strel) that is a line 10.5 µm (21 pixels) long (see inset in Fig. 5(b)), in order to extract the fibrils in the Strel direction. Letting the Strel rotate, a collection of openings by this rotatable line results in a stack of images, each corresponding to one Strel direction. The pixel per pixel maximum in the stack then provides an enhanced image of the collagen fibrils in the slide (Fig.  5(b)). In addition, the Strel direction realizing this maximum provides the fibril orientation. We thus obtain for each pixel of the SHG image the local orientation of the collagen fibrils, which is color coded using the same RGBR scale as for Mueller images in Fig. 5(c). It is worth noting that there is a compromise on the length of the segment. A longer segment gives more accurate measure since more angles are available, nonetheless the length is upper-bounded by local curvature of the fibers inside which the segment needs to fit. tal. This visual impression is confirmed by the orientation histograms displayed in panels (e) and (f), which both peak at 60 • . However, the width (FWHM) of the orientation histogram is clearly larger for SHG than for MP (50 • vs 35 • ). These trends are highly reproducible, as shown in Fig. 7, which summarizes the results of the comparison between the two techniques on 12 ROIs from two different slides. Figure 7(a) displays the peak values of the fibers orientations obtained by the two techniques, which show an excellent agreement. Linear fitting provides indeed a slope of 1.01, very close to unity, with a low intercept angle (0.24 • ) and a high R 2 value (0.93). Nevertheless, Fig. 7(b) shows that the distribution width of fiber orientations obtained by SHG is typically twice the one obtained by MP.
This difference in the orientation distribution widths is not surprising because the two techniques are different both in their principles and in their implementations. First, the spatial resolution was poorer for the polarimetric setup than for SHG (2 µm vs 0.4 µm), which may induce a spatial average of the orientations obtained from MP resulting in narrower distributions. Second, and most importantly, the comparatively larger peak in the SHG orientation histogram is explained by the preliminary mean filter. Indeed, a mean filter in the image domain followed by the linear opening used for detection of orientation results in a convolution in the domain of orientation. The width of this convolution is proportional to arc tangent of the ratio of the size of the linear filter window and of the length of the rotating segment. In contrast, the MP images are not spatially filtered.

Conclusion
In this work, we investigated the possibility to use Mueller polarimetric imaging coupled with Picrosirius red staining to characterize the in-plane azimuthal orientation of collagen fibers in histological slides. The optical anisotropy, and more particularly the retardation, was analyzed and the fibers were assumed to be parallel to the slow axis. Second Harmonic Generation imaging was taken as a "gold standard", due to its high specificity for fibrillar collagen, and its capability to directly provide spatial images of the fibers. Twelve different ROIs, taken from two different slides, were imaged by both techniques. The distribution widths obtained by the two techniques are different, with SHG values systematically larger than MP ones. This discrepancy may be due to different spatial resolutions of the two setups, as well as to SHG image processing, including preliminary mean filtering that widens orientation distributions. In con- trast, the peak values of the distributions obtained by the two techniques were found to be very close meaning that the dominant orientation of collagen fibers, if any, can safely be determined by MP imaging. As MP is much cheaper and simpler to implement than SHG, this technique may be the most suitable technique for systematic studies of histological slides of healthy and diseased connective tissues.