Translational label-free nonlinear imaging biomarkers to classify the human corneal microstructure.

Diseases that affect the cornea can lead to severe vision loss and have tremendous social impact. These diseases are associated to deviations from normal structural order and orientation of collagen fibril bundles. Unfortunately, resolving non-invasively the corneal collagen structure is not possible to date. In this work, polarization sensitive second harmonic generation (pSHG) microscopy is used to obtain information with molecular specificity on microstructure of human corneas. This information is used to develop a set of label-free imaging biomarkers that were generated by means of a novel methodology based on mathematical tensorial calculus. The method is proven to be highly sensitive and robust. The use of these biomarkers permits accurate characterization of the anisotropic, depth-dependent, structural organization of corneal collagen fibril bundles without any a priori information. The method can be valuable to improve understanding of microstructural pathophysiological changes of the human cornea close to in vivo conditions.


Introduction
The cornea is the main refractive element of the human eye. Its main structural component is type I collagen. Diseases affecting the cornea represent a major cause of blindness worldwide, second only to cataract in overall importance [1]. Unfortunately, state of the art diagnostic devices, such as confocal reflectance microscopy or optical coherence tomography, are unable to resolve collagen fibril bundles and cannot provide information on their threedimensional (3D) organization.
The shape and transparency of the cornea are maintained by the highly organized assembly of stromal collagen. These are embedded in an optically homogeneous ground substance, the stromal matrix, consisting of glycosaminoglycans and proteoglycans [2]; in addition, they show a highly heterogeneous distribution over the stromal depth. Collagen fibril bundles in the most anterior stroma show pseudo-random spatial distribution; across the posterior stroma, they are arranged in lamellae that run almost parallel to the corneal surface, forming a grid-like lamellar structure with more ordered organization if compared to the anterior stroma. Clinical benefit of resolving the microstructure of the corneal stroma would extend to several applications, such as screening of keratoconus (the first cause of corneal transplantation in young people worldwide), precise understanding of corneal scarring (first cause of corneal blindness worldwide) and monitoring of the effect of corneal surgery (excimer laser surgery, femtosecond laser assisted lamellar corneal transplantation) [3,4].
Most of the work performed to characterize the microstructure of the human cornea has been carried out by using various optical techniques, such as electron microscopy, circular polarization microscopy, full-field optical coherence tomography, small-angle light scattering and x-ray diffraction [5][6][7][8][9]. Multiphoton microscopy techniques, such as two-photon excited (auto) fluorescence (TPEF), and second and third harmonic generation (SHG and THG, respectively), can be an alternative to these traditional methods, since they provide intrinsic tissue optical sectioning capabilities combined with deep penetration and multimodality signals approach [10]. In particular, polarization sensitive SHG (pSHG) microscopy is a unique tool to reveal structure and orientation of several biological tissues at a molecular level [11,12]. Feasibility of in vivo imaging and analysis of collagen tissues based on spatial information retrieved by pSHG images has been demonstrated previously [13][14][15][16][17][18][19]. Data retrieved by corneal pSHG images can be valuable to develop label-free imaging biomarkers for evaluating the complex microstructure of the corneal tissue.
In this work, we demonstrated the use of pSHG microscopy for providing valuable descriptors of the human corneal structure based on the acquired images; these images were processed using methods based on a tensorial mathematical approach taken from liquid crystal theory. As a result, three parameters were defined that described the architecture organization of the collagen molecules. These parameters are the collagen fibril preferred direction, f, the scalar order parameter, S, and the inhomogeneity index, T. Using these optical imaging biomarkers, the structural order and distribution of the corneal collagen fibril bundles was accurately characterized and their 3D microstructure was automatically classified according to the depth, without any a priori information.

Donor eye-bank tissues
Eight human donor corneas, unsuitable for transplantation, were obtained from the Veneto Eye Bank Foundation (Venezia Zelarino, Italy). All samples were used 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 average age of donors was 58.3 ± 12.1 years, the average post-mortem interval was 12.5 ± 6.1 hours and the average endothelial cell density was 2057 ± 341 cells/mm 2 . The tissues were shipped to the laboratory (Super-Resolution Light microscopy and Nanoscopy Facility, SLN, ICFO, Castelldefels, Barcelona, Spain) via express air courier and preserved at 4° C in corneal storage medium enriched with 15% dextran. Since the microscope objective used in this experiment has a working distance of 200 µm, before shipping, each donor cornea was cut, by the eye-bank technician, to a targeted depth of 220 µm from the anterior stromal surface (anterior lenticule) using a microkeratome (8.0 mm diameter, Moria evolution 3, Moria, France) after gently removing the epithelium with a Merocel sponge (Medtronic, Minneapolis, MN, USA). The posterior lenticule of each corneal tissue was then trephined to a size of 8.0 mm diameter with a hand trephine (Micro-Keratron, Geuder AG, Heidelberg, Germany).

pSHG microscopy imaging
The pSHG microscopy setup has been fully described before [15,19,20]. In this experiment, a Kerr lens modelocked Ti:sapphire laser (MIRA 900f, Coherent, France) was used as the excitation source. This pulsed laser source was operated at a central wavelength of 810 nm, pulse width was 160 fs (measured at the sample plane) and repetition rate was 76 MHz. An oil immersion 1.4 NA 60x objective (Plan Apochromat, Nikon Instruments Europe, Netherlands) was used to focus the light on the sample, while a 1.4 NA light condenser was used to collect the signal generated in the forward direction. Typical frame acquisition time for a single 500x500 pixels image was about ~1.5s. The effect of depolarization of the fundamental beam introduced by the different optical components was also measured at the sample plane, as previously described [12].
Each corneal lenticule was placed in a custom-made cell filled with 15% dextran. The sample was mounted between two different #1 cover-slips (0.13 to 0.16 mm thick), with its anterior surface parallel to the scanning plane. Different z-stacks of images were obtained from the corneal lenticules according to the following settings: images separated 2 µm at depths from 0 to 30 µm, images separated 5 µm at depths from 0 to 60 µm and images separated 10 µm at depths 60-150 µm. Imaging was performed on several areas of the central 2 mm of each corneal lenticule. 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 measurement; any alteration would have been observed clearly as a consequence of photodisruption.
Nine different pSHG images for each plane of focus were obtained by exciting the corneal lenticule with the same number of different linear polarizations [12,17]. Polarisation of the excitation beam reaching the sample was rotated by means of the half wave plate. SHG images were obtained for polarisation values ranging from 0 to 180° at 20° steps. These images were stored for post-processing. The maximum and minimum transmitted pSHG signals were detected when the incoming linear polarization were parallel and perpendicular to the collagen fibril axis, respectively. In an attempt to reduce noise in the images, each polarisation experiment was repeated four times. The overall acquisition time of a pSHG experiment was 1.5 minutes for each focal plane (37.5 sec for each polarisation session).
The generalized biophysical model assumed in this paper is based on that proposed previously [15,20-24].

Development of label-free imaging biomarkers of corneal collagen fibril bundles
pSHG microscopy was used to determine the in-plane orientation of the collagen fibril bundles that are able to generate SHG signal. Several analytical models have been used in the past to characterize the orientation of collagen fibril bundles using SHG or pSHG images of corneal stroma. These included fast Fourier transform algorithms [25-27], Radon transform of filtered FFT image [14], estimation of the azimuth and elevation angles of fibrils [18] and fitting of polarimetric data based on the tensorial formalism of second order nonlinear optics [17]. In previous studies [28,29] it was attempted to characterize the degree of order of the collagen fibril bundles in tissue by estimating the Fiber Orientational Distribution Function (FODF), and from this data the main direction of orientation was determined. This parameter gives the probability of finding a collagen fibril oriented at a certain angle, φ, contained within the imaging plane, z p . The FODF(φ,z p ) represents the probability density function of finding a collagen fibril bundle at certain angle in the image plane: where the angle φ ranges from 0 to π. For processing purposes, and taking into account that we are dealing with discrete quantities, i.e., pixels in an image, the integral of Eq. (1) can be substituted with summation and FDOF is obtained simply by dividing the number of pixels, which represent the collagen fibril bundles possessing a specific orientation, by the total number of pixels representing all the collagen fibril bundles' orientation in the focal plane. Previous approaches based on the probability density function aimed to define the main orientation by calculating the mean, mode or centroid of the distribution of collagen fibril bundles angles; they furthermore aimed to evaluate only the global degree of order by calculating the standard deviation of the mean of FODF or the angular length around the main orientation that includes 50% of the total number of collagen fiber bundles [28,29]. These methods could provide reliable information on the collagen fibril bundles organization only in the case they are highly oriented along the main orientation (i.e., the FODF shows a narrow peak). On the other hand, these methods are inaccurate and provide unreliable results in the case of low oriented networks of fibril bundles, which yield a broader FODF with multiple peaks. This occurs frequently across the human corneal stroma, as it is described in the literature [2,5,6,[30][31][32]. In addition, the FODF fails to provide any local spatial information about each collagen fibril bundle; it provides only a global evaluation of the degree of order around a certain direction but is unable to measure to which extent each fibril bundle is ordered along the main direction. In this paper, a novel approach to solve this problem is disclosed. On the scale of SHG imaging resolution, we consider that the cornea can be represented as a system formed by elongated microstructures with cylindrical symmetry, which are the collagen fibril bundles. These microstructures present a pseudo-random order alignment and can be investigated by the theory developed to describe long-range orientational ordered system, such as liquid crystals [33-35]. If the notation of liquid crystals is extended to corneal collagen bundles, the structural properties of such anisotropic biological system could be described by the system director, f, which is the main direction along which the microstructures are aligned and by the scalar order parameter, S, indicating the extent to which these microstructures are aligned with respect to f [35,36]. The parameters f and S can be used to describe the global measures of orientation and structural order of the collagen fibril that lie within a certain plane z p .
Since the distribution of collagen fibril bundles depends on their depth within the tissue, it is expected that the parameters, f and S, may also vary as a function of z p . According to the liquid crystal formalism, if the orientation of each collagen bundle is known, f and S can be retrieved from a second rank order tensor [36,37] defined as follows: where ⊗ i i u u is the dyadic product of vectors u i , which represent the axis of orientation of the i th collagen fibril bundle lying on z p , Nm is the total number of the fibrils on z p , and I is a diagonal matrix. The Q tensor order approach does not imply any a priori assumption on the orientational distribution of collagen fibril bundles. It is based on well-established theoretical framework of ordered system, such as in the liquid crystal phase, and gives both global and local order information. Using this mathematical approach, the main orientation of an ordered system and its global amount degree of orientation order are directly calculated by diagonalazing the matrix Q; moreover, the local spatial degree of order information of each collagen fiber bundles is retrieved by calculating the second order of Legendre polynomial. The adopted definition of Q(z p ) is coherent with the Saupe order theory used to describe the long-range orientational order of liquid crystal system [35]; Q(z p ) is a traceless and symmetric matrix and thus it can be diagonalized. Therefore, the collagen fibril bundles organization can be described by considering its eigenvectors and associated real eigenvalues. Since the collagen fibril bundles have cylindrical symmetry, the candidate for the preferred direction of the system is the eigenvector associated with the dominant eigenvalue (i.e., the eigenvalue of the greatest absolute value). We term that dominant eigenvalue S, and its associated eigenvector f (that represents the main direction of the collagen fiber bundles). The dominant eigenvalue S ranges from 1 (for a system of high ordered fibril bundles oriented along f) to −0.5 (for a system of fibrils oriented perfectly perpendicular to the director f). When S is zero, the collagen fibril bundles are randomly oriented, i.e., they do not show any preferred direction within the focal plane z p ,in this case the whole image.
The local spatial degree of order of each collagen fibril bundle with respect to f can be also calculated by using the second order of Legendre polynomial: where the brackets denote an average around the first neighboring of the collagen fibril bundle u i , f•u i is the scalar product of the vectors f and u i , the latter defining the axis of orientation of the its single collagen fibril bundle lying on z p . Similar to S, the scalar product of f and u i , s, also takes values from 1 to −0.5 and can be used to graphically map the spatial order of collagen fibril bundles on a local scale. In this model, if the orientation of the collagen bundles in the corneal stroma is known, their distribution in the tissue can be quantified over the 3D volume that they occupy, and hence their structural order can be characterized. Thus, to determine, the x-y spatial angular orientation of the collagen fibril bundles u i (cos(φ), sin(φ)), lying on the focal plane, z p , pSHG images at that plane were acquired from the anterior and posterior corneal lenticules of each sample and processed. This procedure resulted in the determination of the in-plane angle of orientation of the collagen bundles, φ, for every pixel on the focal plane images. For each of these planes, f and S were determined by solving Eq. (2); then, for each pixel on the plane, s was determined by using the second order of Legendre polynomial. An additional structural parameter, the inhomogeneity index, T, representing the percentage of collagen fibril bundles with a negative s value, was introduced. The higher the T index, the lower the structural order of the corneal collagen fibril bundles lying within each focal plane. In addition, polarization sensitive SHG microscopy can be also used to determine the orientation of the pSHG active molecules that are not contained on the imaging plane, but they form an angle with that plane. In order to investigate the off-plane angle of the collagen fibril bundles, δ, the anisotropy parameter, ρ, was extracted from the pSHG images by processing the data.
We assumed that, in the laboratory the coordinate system is X-Y-Z, the laser is propagating along the Z-axis and its linear polarization can be rotated in the X-Y plane in an angle α with respect the Y-axis. The coordinate system of the collagen fibril is chosen in such a way that the y-axis is contained in the X-Y plane and the z-axis coincides with the principal symmetry axis of the collagen fibril ( Fig. 1). With this geometry, the relation between the laboratory and the collagen fibril frames of reference is given by two Eulerian angles, δ (off-plane angle or elevation/transverse angle) and φ (in-plane orientation angle). The detected intensity can be found as [12,17]: where the term ρ is the anisotropy parameter which reflects the polarization SHG non linear response of the collagen fibrils. In particular ρ can be expressed as [17]: 2 2 33 15 cos 3sin where the anisotropy parameter depends on two factors: the first is the effective SHG angle, χ 33 /χ 15 , and the latter is the off-plane angle, δ. For processing purpose, Eq. (3) can be restated as: As inferred from Eq. (4), variation of ρ reflects the anisotropy SHG response of the collagen fibril bundle due its spatial micro-structural heterogeneity. Furthermore, this variation can also be related to the induced distortion of the polarization signal due to the interweaving of collagen fibril bundles running at different orientations and to the tilted off-plane angle δ. Two consecutive processing steps were performed in order to quantify the two factors in Eq. (4). Firstly, we assumed that δ = 0 (i.e., the collagen fibril bundle axis is parallel to the focal plane), so that Eq. (4) reduces in ρ = χ 33 /χ 15 . This means that any spatial variations of ρ are attributed to the molecular changes of the SHG helical pitch angle θ e through: where θ e can be interpreted as the SHG angle of β ννν associated to the collagen's molecule with respect to the symmetry axis of collagen fibril, which coincides with the collagen's triple helix angle. The ratio χ 33 /χ 15 was retrieved by a log normal fit of the distribution of ρ throughout the image. This estimation was then used to calculate the off-plane collagen fibril bundle angle δ through Eq. (4).

Statistical analysis
Three-dimensional stacks of images were generated from anterior and posterior lenticules of human corneas by means of pSHG microscopy. From these images, the orientation of collagen fibril bundles, φ, were determined for each pixel in the stack. Using this only parameter, three label-free imaging biomarkers were calculated according to what has been already detailed in the previous section. These biomarkers included the main fibril director f, the order parameter S and the inhomogeneity index T. Statistics were performed using the SPSS software (SPSS Inc., version 17.0) by one of the authors (ML), who was masked to the corneal lenticules. The normal data distribution of those parameters was verified using the P-P plot within the software. The analysis of variance was used to test significance between S and T values of the anterior and posterior corneal lenticules in each tissue. The Pearson coefficient was used to estimate linear correlation between these variables. Logistic regression analysis [38] was performed, using the parameters S and T as explanatory variables, to classify the structural order of collagen fibril bundles across the anterior and posterior corneal lenticules. Statistical significance was set at P<0.05 for all the tests performed. The main steps of the processing procedure are summarized in Fig. 2.
Computations did not require specific hardware. Fig. 2. Workflow of the calculation done to retrieve the orientation and structural order of the collagen fibril bundles in pSHG images of the human cornea and to predict their membership in the anterior or posterior corneal stroma.

Method robustness
There are several sources of noise that could influence the experimental results, including polarization distortion due to axial polarisation components introduced by high numerical aperture objectives, birefringence of the cornea, diattenuation of SHG signal, cornea absorption, scattering, dispersion and optical aberrations [11,[39][40][41][42]. It was object of this work to quantify the influence of the above sources of noise in pSHG images of the corneal stroma before further analysis. For each polarisation of the excitation beam, the signal measured at each pixel from the acquired pSHG image was considered as the sum of the deterministic generated SHG signal, I SHG (x,y,z), and the noise. Scanning systems are characterized by reduced Speckle noise, and phase matching effects could be ignored here, since the focused spot is smaller than the wavelength of light used. Therefore, it seems safe to assume that noise is a stationary stochastic process, and the image signal can then be expressed as: First, a region of interest (ROI) was selected in each pSHG image of the stack in order to ensure that the features selected were only due to background noise. The ROI was in average 5% of the whole pSHG image. The process N(x,y,z) was then estimated by evaluating the mean and standard deviation of the intensity levels of the ROI and was used to define the spatial noise with its Probability Density Function (PDF). By estimating the PDF, it was found that the noise could be approximated by a Gaussian function with zero mean and variance σ 2 , meaning that N(x,y,z) could be considered white noise with uniform power spectral density equal to σ 2 [34]. The signal power can be expressed as: where c i is the i th -spectral Fourier transform component, whose index i ranges from 0 to 8 since 9 different polarisations were measured. Note at this point that, due to the nature of The intrinsically interwoven nature of the collagen fibril bundles, particularly in the anterior portion of the stroma, may introduce polarization distortion, which ultimately can influence the exact assessment of the orientation of crossed collagen fibril bundles. A threshold based on the SNR of each pixel was therefore used to filter meaningful data. This allowed accurate determination of the preferred orientation of collagen fibril bundles within each focal plane and the anisotropy parameter. Attenuation and scattering were neglected in this experiment, since they have been shown not to contribute significantly to decrease corneal light transmission at 810 nm wavelength [43,44]. The signal intensity was found to exponentially decrease with increasing depth with a characteristic length of ~300 µm; this value was greater than the thickness of each corneal lenticule. In this work, the signal intensity loss was compensated by increasing the excitation laser when possible to generate images with overall constant intensity at different depths. The SNR signal retrieved by experimental data was in the 2 -10 dB range and decreased with depth. Figure 3 shows some results on the distribution of orientations of the collagen fibril bundles depending on the threshold used for the SNR.
The results do not change with the range of threshold values used for the SNR and thus indicate that the range of threshold values used is adequate for these experiments. The fibril orientation varies from 0 (fiber aligned along the horizontal axis, blue colour), to π (red colour). The peak value of FODF does not change, showing that pSHG data analysis is highly robust with respect to the noise. The fibril direction is referenced to the horizontal axis in counter clockwise direction.

Results and discussion
The highly anisotropic, depth-dependent, 3D microstructure of the human cornea is represented in Fig. 4. The novel pSHG imaging biomarkers, f, S, and T were calculated for each plane within each corneal sample. In addition, the FODF was calculated and used as control with the novel methodology to describe the stromal microstructure. The FODF distribution was obtained from φ through pSHG data analysis ( Fig. 5(a)). However, the mode of the FODF did not provide any significant information about the complex spatial arrangement of collagen fibril bundles in the samples (Fig. 5(b)). On the contrary, the novel pSHG imaging biomarkers provided relevant information on the structural order and orientation of collagen fibril bundles according to their depth (Fig. 5(c) and 5(d)).  It provides reliable measurements when only one or two preferred directions run across the focal plane, as found in the posterior stroma; nevertheless, it does not provide any local spatial information to classify the complex microstructure of the human cornea with respect to depth. (c) The main director f of corneal collagen fibers is calculated by taking into account the direction of the collagen fibrils (u i ) that lie within the focal plane z p ; here, f is 101° and 120° for the anterior and posterior pSHG images respectively. Colour scale bars show radians. The fibril direction is referenced to the horizontal axis in counter clockwise direction. (d) The spatial structural order of corneal collagen fibrils is defined by S and can be spatially mapped by s, that ranges between −0.5 (blue colour) and 1 (red colour). The maximum positive value indicates that the collagen fibril bundles are aligned along f; the maximum negative value indicates that they are oriented perpendicular to it. Here, the parameter S is 0.32 and 0.80 for the anterior and posterior pSHG image respectively. The parameter T, indicating the % of collagen fibril bundles with a negative s value within the focal plane, is 29.5% and 0.6% for the anterior and posterior pSHG images respectively.
The calculated value of the global structural order value, S, was overall lower for the anterior stroma compared to that for the posterior stroma in each tissue ( Fig. 6(a)). The average S value of the anterior corneal lenticules were 0.49 ± 0.03 (95% Confidence Interval: 0.48-0.50), while in the posterior ones it was 0.57 ± 0.04 (95% CI: 0.55-0.59); the difference between the anterior and posterior S values was highly statistically significant (P = 0.002). The average T values of the anterior and posterior corneal lenticules were 21.0 ± 0.8% (95% CI: 20.2-21.8%) and 16.3 ± 1.1% (95% CI: 14.9-17.2%) respectively. The average T value was higher for the anterior stroma than for the posterior stroma in each tissue (P = 0.004; Fig.  6(b)). These results can be directly related to the spatial distribution differences between anterior and posterior stroma already detailed in the relevant literature. The anterior stroma presents a pseudo-random distribution of collagen fibril bundles, which translates in a lower S value and a higher T value compared to the posterior stroma, where collagen fibril bundles tend to arrange in two main perpendicular directions. However, each parameter by itself was not enough to discriminate entirely between the anterior and posterior corneal stroma microstructure. In particular, S values on the posterior stroma for samples 3 (0.53) and 5 (0.52) were slightly lower than that for the anterior stroma in sample 2 (0.54); values of T on the anterior cornea for samples 2 (17.5%) and 7 (18.9%) were lower than that for the posterior cornea on samples 3 (19.3%) and 5 (20.5%). Since the S and T variables were highly correlated both in the anterior (r = −0.88, P = 0.004) and posterior (r = −0.99, P = 0.0004) corneal stroma, a logistic regression analysis was performed to define statistically the specific structural organization of the anterior and posterior corneal lenticules of each tissue using both these parameters as descriptors. Cox and Snell R 2 index, from regression logistic analysis, gained a maximum value of 0.75, with Nagelkerke R 2 = 1.0. The statistical analysis thus showed 100% accuracy to discriminate between the anterior and posterior stromal collagen fibril bundles micro-structural organization in each tissue (Table 1). These results clearly demonstrate our ability to quantify the heterogeneity of the collagen fibril bundles distribution and organization of the central portion of the stroma across the corneal depth by calculating the three label-free imaging biomarkers from pSHG images. Biomarker f provided information on the preferred orientation of collagen fibril bundles that lie within each focal plane. The biomarkers S and T provided information on how the collagen fibril bundles are arranged and aligned with respect to the main fibril director f within each focal plane over the corneal depth; S is related to the overall degree of orientation of the collagen fibril bundles in the tissue, while T is related to the percentage of fibril bundles aligned perpendicularly to the main director. The complementary use of the biomarkers S and T as descriptors in a logistic regression analysis of the corneal stroma structural organization was shown to be able to classify unequivocally the degree of order of collagen fibril bundles across the anterior and posterior corneal stroma. The method was used to discriminate automatically, with 100% accuracy, between the fibril organization of the anterior and posterior central stroma. Overall, the analysis of S and T indicated that the distribution of collagen fibril bundles in the anterior cornea is less ordered than in the posterior cornea, in which the collagen fibril bundles lie mostly parallel to the corneal surface and do not show significant vertical branching, which are specific of the most anterior collagen fibril bundles. This was in agreement with previous findings [1,5,16,[30][31][32]45,46] that report short bundles of collagen fibrils of consecutive layers crossing each other at variable angles in the anterior central stroma and wider and longer collagen fibrils forming lamellae that mostly run along two preferential perpendicular directions in the posterior central stroma.
The 3D organization of the corneal collagen fibril bundles was investigated through the analysis of the off-plane angle of collagen fibril bundles, δ. This was on average 9.5 ± 1.7° (range 6.6-11.4°) and 7.8 ± 1.6° (range 5.6-10.1°) in the anterior and posterior corneal lenticules respectively. Overall, these values are in agreement with the description reported on the 3D organization of the collagen fibril bundles on the cornea depending on depth. Previously [31,32], the transverse angle was found to range between 19° for the collagen lamellae lying immediately below the Bowman's layer to 11° and 8° for the collagen lamellae of the anterior and posterior stroma respectively.

Conclusions
Polarization sensitive SHG microscopy is a label-free optical imaging method that can reveal the complex 3D organization of corneal collagen fibril bundles with micrometric accuracy. All methods reported previously can provide information on the global orientation of the individual corneal collagen fibril bundles. However, they fail to provide measurable structural properties and to describe the depth-dependent corneal anisotropy based on data provided by pSHG images.
In this work, we developed and validated a method to quantify automatically the 3D depth-dependent structural organization of collagen fibril bundles in pSHG images of the human cornea without any a priori information. Using this mathematical approach, both the main orientation of the collagen fibril bundles and their global amount degree of orientation order within the focal plane were directly calculated. In addition, the present method provided information on the local spatial degree of order information of each collagen fibril bundle. The automated analysis of nonlinear pSHG signals was shown to be a specific and sensitive probe of the structural organization of corneal collagen fibril bundles. This provided quantitative information never achieved before in studies of the human cornea and due to its non-invasive nature, the technique represents a significant improvement towards translational research in the ophthalmic and applied biophotonics communities. Detecting deviations from the normal structural order of the anterior and posterior stroma of the human cornea at molecular level may benefit studies on the pathophysiology of the most common corneal diseases and enhance the research of methods to treat these diseases in advance with respect to current therapies. Future work will aim to identify associations between the novel parameters and corneal diseases. Deviations from normal structural order of the collagen fibers may be associated to diseases affecting the anterior (keratoconus) or posterior (Fuch's dystrophy) cornea. Furthermore, in an attempt to avoid tissue dissection, the authors are considering the use of long working distance and high NA objectives. The intrinsically noninvasive nature of this technique enhances its potential to transition from fundamental research to clinical studies. The novel imaging biomarkers hold not only great potential for widening the understanding of healthy structural organization, but also for evaluating pathophysiological changes of the corneal microstructure close to physiological conditions. For example, this methodological approach would be applied to understand and quantify if there are changes of arrangement and/or orientation of collagen fibril bundles in keratoconus [47].
Before this work can be applied to studies in human subjects, it is desirable that the technique is optimized. In particular, there are a number of issues that should be addressed in the future. One of these issues is the time required to generate the pSHG images. The authors have already presented some work in this direction, in a method where the pSHG information is generated on a single scan by means of illuminating the sample using circular polarization, and detecting SHG signal generated with three different polarisations simultaneously on three different detectors. Using this method, acquisition time can be reduced from 1.5 minutes to a matter of seconds [48]. Another important issue to consider in the future is the fact that, in an in vivo application, the pSHG experiments should be performed in the backward direction, since it is not feasible to measure the forward SHG signal in a real eye. In this case, the back scattered SHG has been reported to be of around 10% of the forward SHG signal [49]. This will reduce the SNR in the measurements. However, using the fast method mentioned above, the SNR could be increased by repeating the experiment to average noise, and still reduce the total acquisition time. Finally, safety should be a concern in any in vivo experiment, and should be tested in detail in the future. In particular, there are mainly two questions that are critical in terms of laser damage: retinal and corneal damage. There are existing commercial ophthalmic applications of femtosecond pulsed laser sources like the one we used in these experiments. In particular, these lasers are used in cataract surgery and in corneal flap creation for refractive laser surgery [50]. In both of these applications, ocular tissues are photodisrupted and there is no concern with retinal damage. It is safe to accept that, since the experiments presented here did not generate any changes in the corneal tissues, the laser intensities required are lower than those used in femtosecond laser assisted cataract surgery or corneal flap creation. This is therefore an indication that these experiments could be performed with no retinal damage. However, further studies should be required in the future to ensure the safety of the technique.
In addition, the method proposed can be also applied to study the arrangements of elongated fibrous microstructures and their structural order in all human fibrous tissues, benefiting to applied research in the field of oncology, neuroscience and dermatology [14, [51][52][53].