Identifying crossing collagen fibers in human corneal tissues using pSHG images.

Polarization sensitive second harmonic generation (pSHG) microscopy has been used previously to characterize the structure of collagen fibers in corneal samples. Due to the typical organization of the corneal stroma, the information that pSHG provides may be misleading in points where two different collagen fiber bundles orient along different direction crossings. Here, a simulation that illustrates the problem is presented, along with a novel method that is capable of identifying these crossing points. These results can be used to improve the evaluation of corneal collagen structure, and it has been applied to analyze pSHG data acquired from healthy and keratoconic human corneal samples.


Introduction
The shape and transparency of the optical media of the eye, like the cornea or the crystalline lens, play an important role in the quality of image formation on the retina. The microscopic structure of these biological tissues is closely related to their optomechanical characteristics. In the case of cornea, its shape and transparency are provided by this structure. The corneal stroma is formed by a number of stacked lamellae of type I collagen fibrils, although type VI collagen and proteoglycans can also be found in its structure [1][2][3][4]. Any alteration of this structure, caused by either trauma or disease, can lead to changes in the physical and optical characteristics of the cornea, and result in vision loss [2].
Different techniques have been used to characterize the microstructure of the human cornea in the past. X-ray microscopy and scanning electron microscopy have been proposed to precisely characterize molecular structure of the corneal collagen, such as the calculation of its helical pitch angle [2,4,5]. However, these techniques require sample preparation protocols that are not compatible with in vivo measurements. Recently, multiphoton (MP) microscopy has been used to study the microscopic structure of different types of biological tissue, including human corneal samples. MP microscopy has a great potential for in vivo studies due to the following three main reasons: firstly, different molecules present in a biological tissue can be observed by means of MP microscopy without the requirement of any exogenous contrast agents [6]. Secondly, MP signal is usually generated by means of ultrashort laser pulses at intensity values that are safe for in vivo experiments. Furthermore, the applied excitation wavelength is usually close to the infrared, and exhibits high tissue penetration. Thirdly, MP microscopy features intrinsic axial depth discrimination and reduced photodamage.
From all of the different MP microscopy techniques, second harmonic generation (SHG) microscopy has become very popular for studies of the cornea. This is due to the fact that corneal collagen is packed in the tissue in such a way that it behaves like an SHG active polycrystalline lattice [7][8][9][10]. This characteristic has been used in the past to observe the orientation of the collagen fibers using different theoretical models [9,[11][12][13][14][15]. In particular, the application of polarization sensitive SHG (pSHG) has been reported for the classification of human corneal images according to the depth in the tissue at which they were acquired [16]. In these experiments, the polarization direction of the excitation beam is changed, resulting in a modulation in the intensity of the SHG signal obtained from the tissue. This modulation is related to the angle between the polarization of the excitation beam and the orientation of SHG active molecules. However, in this model it is assumed that all of the SHG active molecules within a pixel volume, or voxel, are oriented along the same direction. This requirement can produce misleading results in the study of corneal collagen tissue, since the corneal collagen fibers usually intertwine and are oriented along different directions. This situation is found especially in the anterior portion of the stroma, which plays the most important role in bearing stress and maintaining corneal shape [3,[14][15][16].
In this work, the effects that collagen fiber crossings can have on the results of the pSHG model are illustrated. A method that is aimed at improving the information that can be generated by means of this pSHG biophysical model is also presented. This method is able to detect the pixels of the pSHG images that include the contribution of bundles of collagen fibers oriented along different directions. Once these pixels are detected, they can be filtered out, providing more reliable information related to the orientation and molecular structure of the collagen fiber bundles. These results have been experimentally tested on starch samples and also on healthy and keratoconic corneas.

Biophysical model
The biophysical model used in this work has been described extensively in previous studies [17][18][19][20][21][22]. This model can predict the intensity of SHG signal, I SHG , generated by an active supramolecular assembly with cylindrical symmetry depending on the excitation beam polarization orientation, α, and also the orientation of the SHG active molecule, φ, as follows [17]: where a 0 , a 2 and a 4 are coefficients that will be described in more detail in the following paragraphs. By taking the Fourier Transform (FT) of Eq. (1) with respect to α the following expression is reached [11]: where Ω is the spatial frequency in the Fourier domain, and c.c. indicates the complex conjugate. From Eq. (2), it is possible to calculate two different quantities related to the molecular structure of the SHG active assembly: the orientation of the supramolecular assembly, φ, and the orientation of the hyperpolarizability tensor dominant axis, θ e . In the case of corneal collagen, φ is usually related to the orientation of the collagen fiber bundles, and θ e has been previously related to the helical pitch angle of the collagen triple helix. These identifications will be assumed from now on in this text. However, the general meaning of φ and θ e should be considered when extending this model to other types of tissue.
The helical pitch angle, θ e , can be calculated from the parameters in Eq. (2) as follows [11]: In the case of the orientation of the collagen fibers, φ, it can be determined by computing the complex argument of the second coefficient in Eq. (2) It has also been detailed in previous studies that the values of the coefficients a 0 , a 2 and a 4 can be calculated from θ e, [11,23], and their explicit expressions are: Note that A is a proportionality constant that appears in all of the components and here it will be set to 1, since it won't affect the calculations in Eqs. (3), (4) or (5). According to this, for any given value of φ and θ e , the intensity of the SHG signal can be determined for a particular value of α using Eq. (1). However, the described model assumes that the orientation of the SHG active molecules is well defined, and therefore all the molecules are oriented in a particular direction. As already mentioned, this may not be accurate for collagen fibers in different lamellae of the cornea.
Previous reports have shown that the overall SHG intensity at a pixel where two collagen fiber bundles cross can be described as the sum of the signals generated by each of the isolated fibers [14,15]. This can be expressed as: where I TOT (α) is the total SHG signal intensity generated at a particular point in the sample, illuminated by a laser beam polarized along the direction α, while I SHG,1 and I SHG,2 are the SHG intensities generated on the fibers 1 and 2 respectively when considered independently.
To evaluate the error incurred by the pSHG model used when collagen fibers oriented along different directions are considered to be aligned along one specific and well-defined direction, we generated a set of numeric data. Two different matrices with 512x512 elements were generated, each matrix representing the orientation of a collagen fiber, φ, on each pixel. These matrices are shown in Fig. 1. The values of φ 1 , showing the orientations of collagen fiber 1 from the first matrix relative to the horizontal axis, are shown in Fig. 1(a). The value of φ 1 changes between −90° and 90°. The second matrix with the values of φ 2 , is shown in Fig. 1(b). It has been divided into two halves: in the left half, φ 2 was kept constant at a value of 0°, while in the right half, the value of φ 2 was changed as the angle with the horizontal axis, in the same way as φ 1 . The values of φ 1 and φ 2 have been used to calculate I SHG,1 and I SHG,2 respectively, in Eq. (9). For that, we chose a value of θ e = 45° close to the existing data for collagen [13,18,24] and by applying this value in Eqs. (6), (7) and (8) the values for the coefficients a 0 = 2.875, a 2 = 1.5 and a 4 = −0.375 were calculated. Similar to the actual pSHG experiment, which will be described in following sections, nine different values of α (from 0° to 180° with steps of 20°) were used to calculate I SHG,1 and I SHG,2 using Eq. (1). Finally, 9 different matrices of I TOT (α) were calculated, one for each α value, as the sum of I SHG,1 and I SHG,2 , according to Eq. (9). In this way, a simulation was designed to study the behavior of the pSHG model when two crossing collagen bundles are considered. A 3D set of I TOT (α) data was generated with a resulting matrix of 512 × 512 × 9 elements. In this data set, it will be possible to differentiate two different situations: on the right half of the images the overall intensity will be the contribution of two collagen fiber bundles aligned in the same orientation, while on its left half the overall intensity will be the sum of the contributions of two fiber bundles oriented along different directions. Since φ 2 is kept constant on the left half of the image, while φ 1 is varying along the horizontal axis, all the possible angles between the two different sets of collagen fiber bundle orientations can be easily visualized in the result. Once this theoretical data set was generated, the FT of I TOT (α) was taken along the polarization direction, α, and Eqs. (3) and (4) were used to recover φ and θ e [11]. The results of these calculations are shown in Fig. 2, and the recovered values will be referred to as φ res and θ res from this point on. In the case of φ res , shown in Fig. 2(a), these values have been calculated using Eq. (4). As already explained, on the right half of the image, the results correspond to two collagen fiber bundles oriented along the same orientation, and this orientation is correctly recovered in φ res . In the left half of the image, the results correspond to two collagen fibers oriented along different direc image, the tw angle φ 1 is dif similar situati As already m in the right ha bundles with between φ 1 an As explain (5) wrap the values of φ res,2 between ± 45° instead. The values of φ res,2 wrapped between ± 45° will be referred to as φ' res,2 from now on, and they are shown in Fig. 3(c). In order to show the similarity between φ' res,2 and φ res,4 , the difference between these two quantities, which is called Δφ res , was calculated and is shown in Fig. 3(d). It should be noted that the values of these results range from −90° to 90°, since the maximum difference between two quantities comprised between −45° and 45° can have an absolute value of 90°.
The results in Fig. 3(d) are interesting, and to better explain them, we have displayed them in a different way in Fig. 4. Figure 4(a) shows the data of Δφ res shown in Fig. 3(d), and in Fig.  4(b) the histogram Fig. 4(a) is presented. This histogram shows 5 different peaks: one peak around 0°, two peaks close to values of ± 90° and two more close to values of ± 45°. To understand what these peaks are related to, let us concentrate on the histograms of the left and right halves of the image separately.
Let's first consider the data of the right half of the image in Fig. 4(a). In this part of the image, all the molecules are oriented along the same direction, and it is expected that Δφ res would be close to zero. These are the data in Fig. 4(e). In the corresponding histogram, Fig.  4(f), it can be seen that the most repeated values for Δφ res are around 0°, which correspond to the pixels where the values of φ' res,2 and φ res,4 are very similar. However, there are also smaller peaks close to values of ± 90°. In the image, Fig. 4(e), the pixels with values of Δφ res close to ± 90° (these are the pixels in maroon and dark blue) are found along the directions ± 45°. These peaks come from the noise introduced by the calculation itself. To illustrate this, let us consider a point with φ res,4 = 44°, and assume that the noise changes the expected value of φ res,4 by a 5% to a different amount of 46°. After phase wrapping, the new value of φ res,4 would be −44°, which once subtracted from φ' res,2 will produce a value of Δφ res = 88°. As mentioned above, this situation will mainly arise around the values of φ close to the phase wrapping angle limit, which is indeed observed in Fig. 4(e). In Fig. 4(c), the results for the left half of the Δφ res image are shown. This data presents crossing collagen fibers along different directions. The largest values of Δφ res appear when the crossing fibers form angles of ± 90° (orange and blue pixels). In the corresponding histogram, Fig. 4(d), there is also a main peak around 0° as in the previous case. However, there are also peaks in Δφ res histogram around the values ~ ± 45°. These values correspond to the pixels where the two collagen fibers are not aligned and, they form a certain angle instead. In Fig.  4(c) it can be seen that these pixels correspond to those positions where the difference between the input values φ 1 and φ 2 is larger than 45°, which is where blue and orange pixels can be found.
After this consideration, it can be seen that the histogram in Fig. 4(b) shows all of the peaks described above: one large peak around 0°, two peaks at ± 90° and two more at ± 45°. Using the data in the histogram of Fig. 4(b) we can determine the pixels with crossing fibers, which would be those included in the peaks around ± 45°. These pixels have values of Δφ res in the ranges of (−53.11°, −22.92°) and (22.92°, 53.11°). This information can be used to identify and filter out the pixels with crossing fibers in order to obtain the collagen fiber orientations more reliably.
In the following sections, these results will be validated using the data acquired from different biological samples.

pSHG microscopy imaging
The results obtained with the simulation described in the previous section have been tested using the experimental data from real samples. To do this, the pSHG microscopy images have been acquired using a setup that has been fully described previously [13,21,23]. A Kerr lens modelocked Ti:sapphire laser (MIRA 900f, Coherent) was used as the excitation source. This pulsed laser source was operated at a central wavelength of 810 nm with a pulse duration of 160 fs (measured at the sample plane) and a repetition rate of 76 MHz. A water immersion 1.05 NA 25x objective (Plan Apochromat LWD*, Nikon) was applied to focus the light on the sample, while another 1.05 NA 25x water immersion objective (XLPlan N, Olympus) was used to collect the signal generated in the forward direction. The theoretical axial resolution of the system was 1 µm [25]. Typical frame acquisition time for a single 512x512 pixels image was about ~1.5 s. The effect of depolarization of the fundamental beam introduced by the different optical components was also measured at the sample plane, as previously described [13].
In a first approach, a starch granule was studied as a representative case for the pSHG model. It has been previously reported that the amylopectin molecules in a starch granule are SHG active, and that they are radially oriented within the granule, without any crossing between them [26]. The results obtained will be analyzed in detail in the following section.
Ex-vivo samples of healthy and keratoconic human corneas were also analyzed using our system and method. The samples were obtained in compliance with the guidelines of the Declaration of Helsinki for research involving the use of human tissue and the experimental protocol was approved by the local ethical committee (Clinic Barcelona, Hospital Universitari, Barcelona, Spain). The healthy corneas were provided by the Veneto Eye Bank Foundation (Zelarino Venezia, Italy) with the request of endothelial cell density ≥ 2000 cells/mm 2. They were stored in 15% dextran-enriched corneal medium storage solution [16]. The corneas affected by keratoconus were harvested from patients undergoing penetrating keratoplasty at the Department of Ophthalmology of the University of Florence (Italy). Immediately after surgery, the corneas were immersed in 2.5% glutaraldehyde solution and shipped to the laboratory via express air courier for pSHG imaging.
For imaging, each corneal sample was placed in a custom-made chamber filled with 15% dextran. The sample was mounted on the microscope between two different #1 cover-slips (0.13 to 0.16 mm thick), with its anterior surface parallel to the scanning plane. Z-stacks of images of the whole depth of the sample were acquired. These images were taken at 5 μm intervals. Imaging was performed on several areas of the central 2 mm of each cornea. The optical power after the objective was measured to be ~15 mW and it was adjusted with imaging depth to compensate for signal attenuation. The power was increased up to 50% in deeper planes in order to cover the full dynamic range of the detector. We made sure that no damage occurred to the tissue during measurements; any alteration would have been observed clearly as a consequence of photodisruption [27].
Nine different pSHG images for each plane of focus were obtained by exciting the corneal lenticule with the same number of different linear polarizations [13,19]. Polarization of the excitation beam reaching the sample was changed using a rotating half wave plate. SHG images were obtained for the polarization values ranging from 0° to 180° at 20° steps. These images were stored for post-processing. In an attempt to reduce the noise in the images, each polarization experiment was repeated four times and then averaged. The overall acquisition time of a pSHG experiment was 1.5 minutes for each focal plane.

Results
In this section, the results obtained from the images of real samples are discussed. The section is divided into two parts: the analysis of the images of starch granules, and the same analysis for the excised corneal samples. Figure 5 shows some representative images of the pSHG images obtained from starch granules. Figure 5(a) shows the average of the 9 SHG intensity images with rotated polarizations acquired as described in the previous section. Figure 5(b) shows φ res,2 , the orientation of the pSHG active molecules in each pixel of the image, obtained using Eq. (4). A radial distribution is observed, as expected [26]. To test the validity of the model, Δφ res was also calculated for the sample. The pixels with Δφ res values laying between (−53.11°, −22.92°) and (22.92°, 53.11°) were identified. This information was used to generate a mask to reject these pixels from the data set, since the simulations show that these are the pixels with SHG active molecules oriented along different directions. The resulting image after these pixels were rejected is shown in Fig. 5(c). By comparing this image with the image in Fig. 5(b), it can be seen that only a few pixels of the data set are rejected. These rejected pixels correspond to molecules oriented close to φ = ± 45°. As explained before, these are special values of φ because they are the limit angles for the phase wrapping of φ res,4 values, where noise plays an important role. These results indicate that, except for the particular orientations φ = ± 45°, our method is not affected by the orientation of the SHG active molecules. Figure 5(d) shows the histogram of θ e values. These data have been fitted to a lognormal function. In the plot, θ e values have been normalized to the peak value of the fitted curve. The mode of the fitted lognormal function, θ e = 36.7°, is in a good agreement with the values of the helical pitch angle for starch granules reported in previous works [24]. Since in this case very few points are rejected by the proposed analysis, these results for θ e do not change considerably when using the described mask.

pSHG imaging of human corneal samples
As already mentioned, the use of the method described here has been extended to the study of the collagen fiber structure of healthy and keratoconic excised human corneal samples. Figure  6 shows some representative results, obtained from an excised healthy human cornea at a depth of 220 µm. Figure 6(a) shows the pSHG intensity image of the stroma calculated as the average of the pSHG images acquired for the 9 different polarization directions of the excitation beam, as described in section 2.2. In this image, the structure of the collagen fibers oriented along different directions can be clearly observed. Figure 6(b) shows the results of φ res,2 obtained for this particular set of images using the pSHG model described in the previous studies. In these calculations, noisy pixels have been filtered out in a first approach by considering appropriate values for signal to noise ratio (SNR) of the SHG signal, and appear in black [16]. In the image, the collagen fibers appear in different colors, depending on their orientation calculated using the method described. The main components of the orientation of the fibers appear in the image in orange-yellow, which corresponds to the orientation angles of 25° to 35°. In addition, some dark blue fibers can be identified, corresponding to the angular range of −60° to −80°. These two main directional bands seem to be clearly represented in Fig. 6(a). However, the white circle in Fig. 6(b) shows some fibers in red, which correspond to the orientation of ~70°. Similarly, the red circle shows some fibers in green, which correspond to the orientation of ~0°. This image illustrates the problem that we discussed earlier, since these orientations do not correspond to those observed in Fig. 6(a). As in the from Fig. 6. procedure foll before and af shown in Fig. the data in ord these two his 6. pSHG data acqui pSHG images a lated using the pS do not correspond , a fiber in red (7 in green (0°) seem s where crossing fi c) shows the filtered out u ed to identify p ows that the red set after apply ent with those s compared ou of the SHG a he pixels of the in a single col sing our new m at have been id Fig. 7(c) corres 7. Comparison of t e collagen fiber us information. c) Pix as a way to determ case of starch These data w lowed with φ d fter applying th 8(c). As in the der to estimate stograms is no ired from a norma acquired. Scale ba SHG model, Eq. (4 with apparent orie 0°) appears to be ms aligned along − ibers have been ide orientation da sing Δφ res info pixels with cros d and green co ying our metho shown in Fig. 6 ur results wit anisotropy para e healthy corne llagen fiber usi method based dentified only spond to only t he pixels where p sing a) anisotropy xels that have been mine the difference granules, the v were also mask detailed before he mask. The h e case of the st e the peak valu ot the same, s al cornea at a depth ar is 10 µm. b) 4). The circle high entations of the fib aligned at around −70°. c) The same entified, including ta after the pi ormation. As e ssing fibers and ollagen fibers in d. However, th 6(a). th those gene ameter propose ea data set show ing the SHG a on Δφ res measu by one of thẽ 4% of the pixe pSHG signal has b parameter inform n identified by onl s between the two values of θ e hav ked using the ).    Data acqu described in t shows an aver of 45 μm can results of φ res shows an area + 70° approxi oriented along fibers inside t angles from t o ~ + 75°. In Fig. 8(a) were of the two hist nction, i.e. the vided by the m . 8(c). The dash in Fig. 8(a) an considering all els with crossin ly measured re rved in the mo simulation deta s the overestim 8. Helical pitch ang ejecting pixels wi bution and the resu uired from ke this article.  8(b)). These elical pitch ang ted curves. Th n 2.1, where o lue ( Fig. 2(b)). bers in normal corn s. Scale bar in a) a lognormal funct rneas have al ome representa where the char wn is similar to G model descr and red colored G intensity im ly −40°. Anoth these fibers in orientations in t e crossing col ed in Fig. 8(b) normalized to xels for each θ e d lognormal fu cture show the s of these logn is reduced to values are com gle for collagen his result is co one of the effec nea, a) before and ) is 10 μm. c) H tion for the data in lso been stud ative data for t racteristic colla o that in Fig. 6 ribed. The wh d areas, corresp mage, Fig. 9(a), her interesting ndicate that the the intensity pi lagen fibers h translates on a hardware configuration change to increase the numerical aperture of the objective used. This can be a limit to the performance of the system. Nevertheless, even with optimized axial resolution, misalignments between the image plane and the orientation of the lamellae in the cornea can occur, and this can lead to areas in the image including intersection between different lamellae and therefore the problem may still arise. In this sense, it seems appropriate to have a method that can determine the pixels that may detect conflicting results in an automated way.
We have therefore presented and tested a new method to determine the pixels in pSHG images which may contain information generated at different crossing collagen fibers, and the obtained results are in good agreement with similar methods available in the literature. The presented model has great potential in the study of the collagen structure in the human cornea, especially in its characterization to determine the differences between healthy corneas and those affected by disease known to have an impact on the collagen distribution, such as in the case of keratoconus. Future studies will aim at developing pSHG-based imaging biomarkers for identifying the number and density of crossing collagen fiber bundles in order to assess the structural integrity of cornea tissue.