Quantization of collagen organization in the stroma with a new order coefficient

Many optical and biomechanical properties of the cornea, specifically the transparency of the stroma and its stiffness, can be traced to the degree of order and direction of the constituent collagen fibers. To measure the degree of order inside the cornea, a new metric, the order coefficient, was introduced to quantify the organization of the collagen fibers from images of the stroma produced with a custom-developed second harmonic generation microscope. The order coefficient method gave a quantitative assessment of the differences in stromal collagen arrangement across the cornea depths and between untreated stroma and cross-linked stroma. © 2017 Optical Society of America OCIS codes: (180.4315) Nonlinear microscopy; (170.6900) Three-dimensional microscopy; (170.6935) Tissue characterization; (170.3880) Medical and biological imaging; (100.2960) Image analysis. References and links 1. Y. Komai and T. Ushiki, “The three-dimensional organization of collagen fibrils in the human cornea and sclera,” Invest. Ophthalmol. Vis. Sci. 32(8), 2244–2258 (1991). 2. M. Abahussin, S. Hayes, N. E. Knox Cartwright, C. S. Kamma-Lorger, Y. Khan, J. Marshall, and K. M. Meek, “3D collagen orientation study of the human cornea using X-ray diffraction and femtosecond laser technology,” Invest. Ophthalmol. Vis. Sci. 50(11), 5159–5164 (2009). 3. L. J. Müller, E. Pels, and G. F. Vrensen, “The specific architecture of the anterior stroma accounts for maintenance of corneal curvature,” Br. J. Ophthalmol. 85(4), 437–443 (2001). 4. A. Caporossi, C. Mazzotta, S. Baiocchi, and T. Caporossi, “Long-term results of riboflavin ultraviolet a corneal collagen cross-linking for keratoconus in Italy: the Siena eye cross study,” Am. J. Ophthalmol. 149(4), 585–593 (2010). 5. N. Morishige, A. J. Wahlert, M. C. Kenney, D. J. Brown, K. Kawamoto, T. Chikama, T. Nishida, and J. V. Jester, “Second-harmonic imaging microscopy of normal human and keratoconus cornea,” Invest. Ophthalmol. Vis. Sci. 48(3), 1087–1094 (2007). 6. G. Wollensak, E. Spoerl, and T. Seiler, “Riboflavin/ultraviolet-a-induced collagen crosslinking for the treatment of keratoconus,” Am. J. Ophthalmol. 135(5), 620–627 (2003). 7. D. Cherfan, E. E. Verter, S. Melki, T. E. Gisel, F. J. Doyle, Jr., G. Scarcelli, S. H. Yun, R. W. Redmond, and I. E. Kochevar, “Collagen cross-linking using rose bengal and green light to increase corneal stiffness,” Invest. Ophthalmol. Vis. Sci. 54(5), 3426–3433 (2013). 8. N. Bekesi, I. E. Kochevar, and S. Marcos, “Corneal Biomechanical Response Following Collagen Cross-Linking With Rose Bengal-Green Light and Riboflavin-UVA,” Invest. Ophthalmol. Vis. Sci. 57(3), 992–1001 (2016). 9. M. Kohlhaas, E. Spoerl, T. Schilde, G. Unger, C. Wittig, and L. E. Pillunat, “Biomechanical evidence of the distribution of cross-links in corneas treated with riboflavin and ultraviolet A light,” J. Cataract Refract. Surg. 32(2), 279–283 (2006). 10. S. Kling, L. Remon, A. Pérez-Escudero, J. Merayo-Lloves, and S. Marcos, “Corneal biomechanical changes after collagen cross-linking from porcine eye inflation experiments,” Invest. Ophthalmol. Vis. Sci. 51(8), 3961– 3968 (2010). 11. N. Bekesi, P. Gallego-Muñoz, L. Ibarés-Frías, P. Perez-Merino, M. C. Martinez-Garcia, I. E. Kochevar, and S. Marcos, “Biomechanical Changes After In Vivo Collagen Cross-Linking With Rose Bengal-Green Light and Riboflavin-UVA,” Invest. Ophthalmol. Vis. Sci. 58(3), 1612–1620 (2017). 12. G. Wollensak and E. Iomdina, “Biomechanical and histological changes after corneal crosslinking with and without epithelial debridement,” J. Cataract Refract. Surg. 35(3), 540–546 (2009). 13. S. Hayes, C. S. Kamma-Lorger, C. Boote, R. D. Young, A. J. Quantock, A. Rost, Y. Khatib, J. Harris, N. Yagi, N. Terrill, and K. M. Meek, “The effect of riboflavin/UVA collagen cross-linking therapy on the structure and hydrodynamic behaviour of the ungulate and rabbit corneal stroma,” PLoS One 8(1), e52860 (2013). Vol. 9, No. 1 | 1 Jan 2018 | BIOMEDICAL OPTICS EXPRESS 173


Introduction
How collagen fibers are arranged can have a huge effect on the behavior of collagenous tissue.For example, the sclera of the eye is opaque, while the stroma is transparent despite both tissues being primarily composed of type I collagen fibers [1].In human corneas, x-ray scattering measurements revealed that collagen fibers have a preferred orientation in the stroma, with many lamellae oriented either along the superior-inferior axis or along the temporal-nasal axis except at the periphery of the posterior cornea, where they have a radial preferred orientation [2].Measurements in scanning electron microscopy show alternating collagen lamellae orientation along the two parallel axes in the posterior cornea, although the lamellae of the anterior stroma appear to be more disorganized [3].A change in the orientation and structure of the collagen matrix, such as in keratoconus, has been proposed to result in corneal biomechanical weakening [4], which leads to changes in corneal shape and a degradation of the optical properties of the cornea.A competing theory on keratoconus suggests a change in the interweaving of the anterior stroma with the Bowman layer is the underlying cause [5].
A common treatment of keratoconus is cross-linking, where a photosensitizer, typically riboflavin [6], and light irradiation (generally UVA light) is used to stiffen the cornea.More recently other photosensitizers, such as Rose Bengal and excitation with green light, have been proposed [7].Mechanical measurements (in corneal strips as well as in intact corneas and eyes) have shown increased stiffness in the cornea following cross-linking in experimental models both following treatment ex vivo [7][8][9][10] and in vivo [11,12].From a microscopic standpoint, the photosensitizer promotes bonds on the surface of collagen fibrils or within the surrounding proteoglycans [13].Macroscopically, cross-linking improves the topography of keratoconic corneas to approach normal corneas six months after treatment [14].However, improvements in organization have thus far been reported qualitatively, and a need for a more quantitative metric is present.
Collagen matrices are a natural generator of second harmonic generation (SHG) [15] due to the geometry of individual collagen fibrils [16].The stroma of the cornea [4], composed mainly of collagen types I, III and V [17], is an example of such a matrix, thus making SHG an attractive option for non-destructively and non-invasively measuring the structure of the stroma [5] while opening up the option for direct 3D mapping of collagen fibers [18] by high resolution measurement of the collagen lamellae [19].In general, only forward scattered imaging of stromal collagen gives images of individual fibers, but backscatter imaging can still give a sense of lamellar orientation [20][21][22].
Images of collagen layers taken from SHG measurements can give a qualitative feel for how the collagen fibers are organized, but steps have been taken to provide a more quantitative value of this organization.To this end, a variety of different imaging and processing techniques have been employed, such as Radon transformation [23], structure tensor analysis [21], multi-polarization nonlinear imaging [24], fiber angle distribution mapping [25], texture analysis [26], curvelet transformation [27], and Fourier transformation [28][29][30][31][32][33].The Fourier transform can be done on separate image planes as a whole [28,33] or subdivisions of the image [29][30][31], as well as in three dimensional volumes [32].The work presented here continues and improves on the methodology present in ref [28], as this work quantifies the organization of SHG stromal images.
In this study, the organization of collagen fibers inside of porcine corneas was compared under three different treatments; untreated, Riboflavin-UVX crosslinked (UVX), and riboflavin only treated corneas.SHG signal from collagen fibers was collected and images of the fiber network were pre-processed and Fourier transformed.A metric to quantize the degree of organization of collagen fibers, based on the direction of the Fourier transform, was introduced.

Samples
A total of 25 porcine eyes were obtained from a local slaughterhouse and used within 48 hours post-mortem.The eyes were divided into three different groups; eyes with the epithelium removed (untreated, n = 8), eyes that were treated with the standard Dresden cross-linking protocol (n = 8) [6], and eyes treated with riboflavin solution (n = 8).The remaining cornea was excised from the globe and placed in 4% paraformaldehyde for four hours, to serve as a control cornea to test the repeatability of the image analysis presented below.Paraformaldehyde was used to preserve the cornea because the repeatability measurements took a substantially longer amount of time to complete (5 hours for repeatability measurements vs 0.75 hours for standard measurements).After treatment (described in the following sections), corneal buttons were excised from the globes and placed on a microscope slide with a drop of buffered phosphate solution (PBS) to prevent sample dehydration, which were then covered with coverslips.A drop of water was then placed on top of the coverslip.

Corneal crosslinking
Eyes were crosslinked using the standard Dresden protocol [6].To summarize, the epithelium was mechanically removed from the eyes and a drop of riboflavin (0.1% riboflavin in 20% dextran solution) was applied to the de-epithelialized cornea every five minutes for 30 minutes.After the initial 30 minutes, the eye was placed under an Ultraviolet-A lamp (0.25 W/cm 2 ) and riboflavin application was continued in the same mode as before for another 30 minutes.Eyes treated only with riboflavin were also de-epithelialized and a drop of riboflavin was applied to the cornea every 5 minutes for 1 hour without any ultraviolet-A irradiation.

Imaging system
Fig. 1.Sketch of two-photon microscope setup.The pulsed source was a MaiTai Ti-Sapphire femtosecond laser.The beam position and angle were controlled by a pair of scanning galvanometer mirrors (SGM), and the beam passed through a pair of lenses (L) before entering the focusing objective (F.Obj.) in order to overfill the entrance pupil of the objective.SHG was collected in the forward direction by a collection objective (C.Obj.) and was separated from the excitation wavelength by a pair of dichroic filters (F).SHG signal was single photon counted with a photomultiplier tube (PMT), and a National Instruments PCI card (NI-PCI).The SGM were also controlled with the same PCI card.A beam block (BB) was employed to absorb excess excitation light filtered from the SHG signal.
A custom built two-photon microscope (Fig. 1) with SHG capability was used to generate second harmonic signal from the collagen fibers in the cornea samples.The microscope used an infrared, pulsed laser source (MaiTai eHP, 800nm, 70fs pulses, 80 MHz, average power of 28.5 mW) for excitation.A wavelength of 800 nm was used since it produces the maximum second harmonic signal for corneal collagen [20].The excitation beam was focused with a 0.9 NA water immersion objective (Olympus LUMPLANFL60XW) and SHG was collected in the forward direction with a second objective (Olympus LUMPLANFL40XW, 0.8NA).The resulting signal was filtered with two dichroic filters (SHG signal passed through a 395-415nm bandpass filter and reflected off of a 665nm long pass filter) and collected with a photomultiplier tube (Hamamatsu H10682-210-PMT).A pair of scanning galvanometer mirrors (Cambridge Tech Model 6215H) rastered the excitation focus across the sample, while a pencil motor (Thorlabs Z825B) moved the sample along the optical axis.

Imaging protocols
A custom Visual Basic code used collected signal from the photomultiplier tube to generate images.Collagen images were collected on 187.5 x 187.5 μm surface with a 500 x 500 pixel resolution at approximately the center of the sample.The dwell time at each pixel was 16 μs, for an approximate collection time of 4 seconds/image, and each pixel had a corresponding scale of 0.375 μm/pixel, which is smaller than the diffraction limit of the microscope objective (0.56 μm).The lateral image was oversampled to prevent collagen fibers from being lost or unmeasured.Initial z-sacks had a step size of 2 μm (n = 15 volumes, 5 volumes in each treatment), but the step size was increased to 2.5 μm (n = 9 volumes, 3 volumes in each treatment) later in order to decrease the total amount of collection time for a stack.The increased step size is slightly larger than the axial resolution (2.4 μm).Z-stacks were taken along the entire depth of the cornea, with depths approximately 1 mm although depths varied between eyes both in different treatments and in the same treatment.A preview image was generated at 100 μm intervals along the optical axis to ensure that images could clearly show collagen fibers throughout the stroma without indications of folding or other external effects introduced from excision or placement on the microscope slide.The beginning and ending point of the z-stack was set 20 μm above the anterior surface and 20 μm below the posterior stroma.Post collection, image stacks were inspected to determine the exact location of the anterior stromal surface and if all images collected were viable for analysis.An image was determined to not be viable if no fibers were visible in the image.All images collected after these images were also unused.

Image analysis
A custom Matlab code was developed to analyze the SHG corneal images.The initial steps of image analysis are in line with the method proposed by Bueno et al. [28], which consist of image pre-processing and Fourier transformation.Besides noise reduction, the pre-processing method also incorporated edge sharpening, which followed the procedure laid out in ref [34].First, the first and last three rows and columns of pixels were removed from the image to prevent any edge effects introduced during image collection from being analyzed.Fiber edges were sharpened by subtracting a low-pass filtered image from the original image multiplied by a constant (experimentally set at 2 for this set of experiments).The low pass filter was performed using a Gaussian mask with square size of 10 pixels and standard deviation of 10.The sharpened images were then Fourier transformed after they were windowed with a 2D Hanning window (window constructed with N-point symmetric Hann window in both coordinates, where N is the pixel side length of the image) to reduce the artificial high frequencies introduced in the transform by finite-length sampling.Then, the 270 elements (points) of the Fourier transform with largest magnitude were saved and the rest of points were discarded.The number of points was kept constant for every image, which allowed direct comparison between different measurements.The number of points was experimentally set at 270 because that number of points was sufficiently large enough that lamellar orientation could be measured.For this specific application, 270 points was sufficient, but for other tissue or image types a different number of analysis points may be more appropriate.Saved points were plotted in polar coordinates, providing magnitude and direction for the fibers.From here, the angles were sorted into a total of N windows = 36 ten-degree windows and the number of fiber orientations falling into each window were counted.Ten-degree windows were selected in this work because windows of this size were narrow enough to contain one lamellar orientation yet wide enough to allow for small deviations in the individual lamellae, present naturally or introduced from relaxation of the cornea post excision from the ocular globe.A visual inspection of the lamellar orientation through representative samples of porcine cornea was undertaken to ensure that neighboring lamellae of different depths did not occupy the same analysis window.Points within 2.8° of the 90° and 270° axis were excluded in order to prevent artifacts introduced during image collection from biasing the analysis.The number of points falling inside each window M i is counted from a total number of M total = 270.The order coefficient was introduced as a metric for quantizing how uniform the direction of collagen fibers ran.The order coefficient is defined as the following: ( ) The order coefficient can be interpreted as the normalized standard deviation of the points across windows.The numerator of ( 1) is a measure of the deviation of the proportion of points falling within the different windows with respect to the mean number of points per window (1/N windows ), and the denominator serves as a normalization factor that ensures the order coefficient has a value between 0 and 1.The Order Coefficient = 1 if all the points fall just in two windows (because of the symmetry of the Fourier transform, this is the lowest number of windows possible), and represents a highly organized and uniform orientation across fibers.The Order Coefficient = 0 if points are equally distributed in all directions (defined windows), and represents a completely random distribution of fibers in the image.Figure 2 illustrates the analysis process from the raw image to the order coefficient calculation.The analysis process with synthetic images is demonstrated in Appendix A. Note that the value of the order coefficient can be affected by the number of points used (more points tend to lower the order coefficient as features apart from the collagen fibers are assigned analysis points) and the size of the angular windows (finer windows tend to lower the order coefficient as fewer points are found within the same window).However, the trend of the peaks and valleys between different layers of the same cornea is preserved but less pronounced when the overall order coefficient is lowered.

Repeatability of order coefficient measurements with a paraformaldehyde cornea
The analysis technique was first demonstrated on a cornea fixed in paraformaldehyde to evaluate multiple measurements on one sample in one place.Paraformaldehyde was used to avoid artifacts due to dehydration, specifically affecting spacing between lamellae.After soaking for four hours, the paraformaldehyde cornea was placed in the microscope and four z-stacks of SHG images were taken with common scanning parameters, i.e. the same starting and ending point along the optical axis with the same step size at the same position on plane.Collagen images were taken as described above in Image Analysis, with a total image depth of 400 μm measured with images collected every 2.5 μm, for a total of 161 images recorded for each volume.A visual inspection of the images in the four z-stacks confirmed that all the images were usable for analysis.The average order coefficient for the four stacks was computed for all images in a stack.

Measurement of the degree of organization throughout the untreated porcine cornea and repeatability of processing
The repeatability of the technique was evaluated by measuring the cornea in two orientations.Untreated porcine corneas (n = 3) with epithelium mechanically removed were imaged in the microscope beginning with the anterior side of the cornea placed closer to the excitation objective.All three corneas were imaged on the same day, which minimized any variation in measurement apart from variation in the corneas themselves, e.g.beam path alignment or longer post-mortem measurement time (all corneas measured <12 hrs post-mortem).Images of collagen were collected as described above.Z-stacks with depths of 1200 μm and step sizes of 2 μm between images in the volume were taken for each of the untreated stromas, producing 600 images in each individual z-stack.After the z-scan finished, the samples were flipped over and posterior-to-anterior z-stacks of images were taken.Care was taken to keep the orientation of the corneas the same, i.e. the superior/inferior axis was parallel to the polarization of the excitation beam.Z-stacks of images were taken with the same measurement parameters used in the anterior-to-posterior z-scans.Order coefficients were computed for each image in the stacks and the average order coefficient across all three samples in 100 µm depth sections was found.

Imaging of cross-linked corneas
The three different treatments, as described earlier in the Samples section, were prepared and imaged in the two-photon microscope as described in the Imaging protocols section, with the exception of one UVX cornea (n = 1) which had an expanded image size of 225 x 225 μm surface with a 600 x 600 pixel resolution and had a step size along the optical axis of 2 μm, with the number of analysis points used per image also set at 270.The remaining UVX corneas (n = 7) were imaged and measured in the normal way.Order coefficients were calculated for the first 700 µm of the stroma or until an unviable image was found in a zstack.The resulting values were then divided into three groups based on depths in the stroma: the first 100 μm of the stroma, the 100 μm to 400 μm range, and the 400 μm to the 700 μm range.Mean order coefficient values and variance of the order coefficient were calculated for each treatment and depth range.

Single sample repeatability results
Figure 3 shows the order coefficient vs depth of a single cornea fixed in a 4% paraformaldehyde solution along the same z-scan axis.The average order coefficient along the same 400 μm length (number of images in each z-stack n = 161) was recorded as 0.269, 0.271, 0.272, and 0.278, with variances of 1.86 × 10 −3 , 1.71 × 10 −3 , 1.82 × 10 −3 , and 1.88 × 10 −3 respectively.The differences between averages in the four separate z-stacks were not found to be statistically significant (ANOVA, p>0.05).As well as similar values of average order coefficient across the scan length, the depth of and relative distance between local order coefficient maxima and minima is common between the separate image stacks, as can be seen in Fig. 3.For example, the maximum peak values between all four measurements were within 10 μm of each other (average value 385 μm below the anterior surface).The global minimum for the first three formaldehyde measurements is located at 250 μm below the anterior surface, with a local minima of the fourth located at the same depth.The distance between the maximum located at 100 μm (in the first three measurements, 82.5 μm in the fourth) and the next maximum is 25 μm further below the anterior surface of the stroma (125 μm in the first three measurements, 107.5 μm in the fourth measurement).

Anterior-to-posterior and posterior-to-anterior stromal measurements
Figure 4 shows the order coefficient averaged over three untreated corneas, with epithelium mechanically removed, and over 100 μm depths (number of images averaged per point in Fig. 4 n = 150).The trend of the order coefficient of the anterior-to-posterior z-stacks is mirrored by the trend of the posterior-to-anterior z-stacks, although the order coefficient of the posterior-to-anterior scans tends to be higher at all depths.The difference in order coefficient in the posterior stroma (700-1200 μm) is statistically significant (ANOVA, p<0.05), while the differences in the anterior stroma (0-600 μm) were found to be statistically insignificant, except for the average of 0-100 μm and the average at 400-500 μm (ANOVA, p<0.05).The slight increase in order coefficient of posterior-to-anterior scans, especially in the posterior cornea, can be attributed to how the excised cornea sits and relaxes on the microscope slide, as the posterior side of the excised cornea has an uneven surface causing a slight deformation in the cornea during anterior-to-posterior scans.Figure 5 shows representative images of the collagen fibers in different depths of the stroma, namely the anterior, intermediate, and posterior stroma.Images of the stroma immediately underneath the anterior surface are characterized by short, quickly terminating fibers, as in Fig. 5(a).In contrast, the anterior stroma 100 μm beneath the anterior stromal surface produced images with longer, thick fibers with lamellae visible, as seen in Fig. 5(b).This is reflected in Fig. 4 with the order coefficient initially lower in the primary 100 μm but climbing afterward.In Fig. 5(c) an image taken from the intermediate porcine stroma is seen, which is characterized by thin collagen fibers and highly interwoven lamellae.This region, in the 500-900 μm range below the anterior surface, has a lower order coefficient than the anterior segment of the stroma.Figure 5(d) shows a representative image of the collagen fibers in the posterior stroma, where the collagen fibers are thicker than in the intermediate stroma and the lamellae are more sheet like, i.e. the interweaving between lamellae is reduced.This type of structure translates into a higher order coefficient, as seen in the 900-1200 μm range of Fig. 4.

Order coefficient value of cross-linked corneas
The differences between the untreated, riboflavin treated, and UVX stroma can be seen visually from the individual images collected.Figure 6 shows images of the different treatment types at a depth of 108 μm below the anterior stromal surface.The untreated and riboflavin only treated corneas had crimps, or waviness, in the fibers that can be attributed to relaxation after excision of the cornea from the rest of the ocular globe.However, the fibers of the UVX corneas were in general straighter after excision due to the stiffening of the fibers from the treatment.Figure 7 shows the order coefficients averaged over 20 μm sections and across all corneas in each separate treatment (untreated n = 8, riboflavin only n = 8, UVX n = 8).The first 100 μm have similar values of order coefficient, as the interweaving of the collagen fibers is especially great in the porcine stroma.After the first 100 μm point, the UVX treated eyes had a rise in order coefficient above the untreated and riboflavin only treated corneas.The order coefficient value of the UVX treated corneas returns to the level of the other treatments at approximately 300 μm below the anterior stromal surface.The different sections of the stroma, the layers here divided into the first 100 μm, the 100-400 μm range, and the 400-700 μm range, have a normal distribution in all treatments and depths, which indicates an oscillation of order coefficient.The order coefficient oscillates because some images capture two or more collagen lamellae, which naturally have a different collagen fiber orientation, in the same image.The analysis points in turn are spread across more windows, lowering the order coefficient.Figure 8 shows the individual distributions for the different depth ranges and treatments, and Fig. 9 shows the individual means from the different depth ranges.The total number of images used in the analysis for each treatment and depth is listed in Table 1.
In the first 100 μm, the mean values of the order coefficients are approximately equal across all treatments, with the riboflavin only treatment having a slightly greater (0.01) but statistically insignificant average order coefficient (difference in mean between untreated and riboflavin statistically not significant, ANOVA p>0.05, between riboflavin and UVX not significant, ANOVA p>0.05, between untreated and UVX not significant, ANOVA p>0.05).
In the 100-400 μm range, the mean value of the UVX treated eyes was higher than in the untreated eye and the riboflavin only eye, which indicates that the images produced from the UVX treated eyes had more collagen fibers aligned in the same direction than the other two treatment types (ANOVA pairwise, p<0.05 for difference of mean between UVX and riboflavin only and UVX and untreated stroma).The riboflavin only treated eyes showed less fiber organization than the untreated eye in this range (ANOVA pairwise, p<0.05).After the 400 μm mark, the average order coefficient value of all three treatments equalizes, with the untreated cornea showing a small but insignificant difference in mean order coefficient (statistical significance between combinations of untreated, riboflavin only, and UVX stromas, ANOVA p>0.05).    ) is listed in Table 1.Black bars represent the variance.

Discussion
The procedure presented here evolved from previous quantification attempts of the organization of the collagen distribution in corneas [28] but takes a step further in quantifying how uniform the directions of the collagen fibers run.This procedure is simple in two ways, in setup and in computation.The setup contains the standard elements of a multiphoton scanning laser microscope.The analysis procedure does not require prior knowledge of the collagen fibers in order to perform the analysis or use of complex formulae in the calculation of the order coefficient.The process is also completely automatic and does not require any input, such as drawing of best-fit lines or selection of regions of interest, from the user.This study also mapped out the entire depth of the stroma, which gives the ability to compare different sections of the stroma against each other without the need to mechanically section the cornea.
Measurements in a paraformaldehyde cornea demonstrated that estimates of the order coefficient are consistent across repeated measurements on the same sample.The weaving between different collagen lamellae in the cornea created the local maxima and minima, which are common across all measurements.The lamellae have an alternating orientation, which can range between 0 and 90 degrees [35], and when more than one lamellae is imaged at the same time, the order coefficient is naturally lower as more than one collagen orientation is imaged, that is, the distribution of analysis points is spread out across more windows than when a single, uniform lamellae is imaged.The location and magnitude of the peaks does differ in some image stacks, likely due to relative shifts between the excitation volume and the sample between image stacks, as the shape of the excised cornea after paraformaldehyde soaking was also preserved, which in turn would lead to shifted image depths.The natural curved shape of the cornea precluded a flat surface for the coverslip to sit on, which possibly caused a shift in the focal volume and/or the cornea.In addition, the precision of the pencil motor was limited to ± 1 μm, which is approximately 40% of the axial resolution of the microscope.Despite these possible sources of error, the order coefficient trend was still consistent and the average order coefficient across all measurements was approximately equal, with no statistically significant difference found between measurements.In addition, the average order coefficient of anterior-to-posterior z-stacks was mirrored by posterior-toanterior z-stacks, further supporting the consistency and repeatability of the order coefficient estimates.
The estimated overall trend of the order coefficient in the untreated porcine cornea is consistent with reports in the literature, showing that the collagen fibers in the anterior of the porcine cornea exhibit more order than the posterior stroma due to an increased packing density with depth [36] and a decrease in interfibrillar spacing in the posterior stroma [22], while the anterior lamellae in the porcine cornea close to the anterior stromal surface are characterized by shorter visible fiber lengths [36].This may be due to transverse collagen fibers in the anterior segment of the stroma, as seen in the human cornea [25], although, to our knowledge, no mention of transverse collagen fibers or lamellae in porcine stroma has been reported.
These random and short collagen fibers in the first 100 μm changed into longer and more ordered collagen fibers deeper (i.e. at depths greater than 100 μm from the anterior surface) into the anterior stroma, as seen in the higher order coefficient average at the 100-200 μm mark in Fig. 4 and when visually comparing Fig. 5(a) and Fig. 5(b), as preferred directions and alternating lamellae began to reveal themselves.The dip of the order coefficient trend can be attributed to the increase of interweaving between collagen lamellae in the intermediate porcine stroma, an example of which is shown in Fig. 5(c).The climb of the order coefficient after the intermediate dip comes from the collagen lamellae becoming less interweaved and more sheet-like, as in Fig. 5(d), which is also seen in human corneas [18], leading to a decrease in the number of lamellae imaged at the same time and thus an increase in measured order coefficient.
Measurements of the order coefficient reveal aspects of cross-linking of the porcine eye model.The first is that the first 100 μm of the porcine eye, which are characterized by short collagen fibers and no preferred fiber orientation, show no change in organization after crosslinking.This depth value may vary across samples, as the start of z-stack analysis depends on where the first complete image of the stromal collagen was taken, and this may be affected by slight tilt or uneven location of the cornea on the coverslip.The effect of UVX was greatest, with respect to the order coefficient, in the 100-400 μm range.In the literature, the depth range of effective UVX depth was in the range of 300-400 μm [9,37], but note that this range includes both the human model and porcine model eye.However, these depth ranges were measured by stress-strain measurements [9], post treatment hydration experiments [38], keratocyte apoptosis [39], and depth of demarcation lines [37].The first three are indirect measurements of UVX treatment, while the relationship between the demarcation line and the microscopic effects of cross-linking on the stromal collagen is not well established.Measurements of the order coefficient confirm that the majority of cross-linking takes place in the anterior 400 μm of the porcine eye by directly observing and measuring the collagen fibers.
Also of note is the lack of interocular pressure for the measured samples.The corneas were excised from the ocular globe in order to facilitate forward scattering second harmonic imaging.This allows the collagen fibers to relax, which produces the waviness in the images collected from the untreated and riboflavin only treated stromas.The waviness of fibers increases the spread of analysis points across different windows when computing the order coefficient, thus lowering the value.However, as the collagen fibers had been artificially stiffened, the UVX treated corneas were able to maintain structure and form post excision without any external force applied to the corneas.The fibers of the UVX treated cornea were straighter and without the waviness of the other treatments, thus producing a higher value of the order coefficient.After the 400 μm point, the UVX images showed relaxation and crimping as well.
From the measurement and calculation of the order coefficient, we quantified the degree of organization, especially linear organization and straightness of collagen fibers in the cornea.The collagen organization in untreated corneas were plotted through the depth of the stroma, with the anterior stroma beginning to transition into the intermediate stroma at approximately 500 μm below the anterior stromal surface and into the posterior stroma at 900 μm below the anterior stromal surface, which mirrors qualitative assessments from previous literature sources.The depth at which UVX has a noticeable effect on collagen structure and organization is also confirmed by measurements of the order coefficient.
In conclusion, we have introduced a novel way of quantizing the organization of collagen fibers in the Second Harmonic Generation forward microscopic corneal images.The technique uses the Fourier transform to map out the orientation of collagen fibers and counts the number of fiber orientations that fall within 10° orientation ranges.This method was used to quantify the arrangement of collagen fibers in an untreated porcine eye and the change in organization of collagen fibers in porcine corneas after cross-linking.With this metric, it was found that the corneas following UVX had a higher degree of organization than both untreated eyes and eyes only treated with riboflavin in the 100-400 μm range below the anterior stromal surface.Future studies can use the order coefficient to evaluate the effectiveness of different cross-linking treatments and evaluate the evolution of collagen organization in keratoconic corneas.

Appendix A. Ideal images and the order coefficient
In Fig. 10(a), the image is Gaussian distributed noise, which represents a perfectly disordered system.The Fourier transform is shown in Fig. 10(d), and has an order coefficient of 0. In contrast, image Fig. 10(b) is composed of alternating white and black bars two pixels wide.The Fourier transform of the bars is shown in Fig. 10(e), which produces an order coefficient of 1. Image Fig. 10(c) combines three alternating black and white patterns in one image, with bar widths of 4, 8, and 16, and the breakdown of analysis point generation was 18%, 26%, and 56% respectively.In general, more analysis points are generated from thicker bars than thinner bars, giving greater weight to thicker bars in Fourier space.The Fourier transform of the image is displayed in Fig. 10(f) and the order coefficient = 0.617, which is higher in comparison to the order coefficient value of an even distribution of points across six windows (0.542).The difference in bar size, as seen from Fig. 10(f) and the resulting order coefficient, produces anisotropy in Fourier space.In terms of imaging in stromal collagen, individual collagen fibers do not tend to vary in size considerably and such variations would be expected to be much smaller than the resolution of a light microscope.Fig. 10(g) is an image of bars with thicknesses ranging from 1-3 pixels wide and spacing 1-3 pixels wide.The same image was rotated 90°, producing equal numbers of points in analysis windows, as shown in image Fig. 10(j), and an order coefficient of 0.686.This value is invariant so long as the analysis points are distributed evenly in four different analysis windows (two windows per direction).Fig. 10(h) shows two crossed 5-pixel patterns with equal intensity, which produces the Fourier transform seen in image Fig. 10(k) and an order coefficient value of 0.686.Fig. 10(i) is two crossed 5-pixel patterns, but one pattern had an order of magnitude less intensity.The Fourier transform, image Fig. 10(l), has more points along the higher intensity (61% in the higher intensity distribution), and the order coefficient had a value of 0.705.This demonstrates that differences in intensity have a small but noticeable effect on the order coefficient by creating anisotropy in the distribution of analysis points.
Fig. 10(m) is an image of alternating white and black bars one pixel wide with the intensity of the white bars randomly varying (2-5 pixel long blocks with values ranging from 0-150 with value steps of 15, poissonian noise added to image) rotated 109°.The corresponding Fourier transform is shown in Fig. 10(p) and the order coefficient was calculated to be 0.388.Fig. 10(n) is another image of alternating white and black bars one pixel wide, but with two perpendicular patterns.The first pattern is the same as in Fig. 10(m) while the perpendicular pattern has section 20-50 pixels long with the same value range and discrete values.The majority of analysis points, as seen in Fig. 10(q), were generated from the pattern with the longer constant pixel sections, and the final order coefficient of the combined image is 0.683, slightly lower than the order coefficient of the longer pixel pattern only (0.707).From the analysis of these two images, it can be seen that patterns of long fibers weigh more heavily in order coefficient analysis than shorter fibers.To review, the order coefficient has the lowest value possible when the distribution of analysis points is evenly distributed and rises when anisotropies are present.The order coefficient, as a scalar function, is rotationally invariant.Values of the order coefficient are affected by bar or feature width in the case that two different widths are present in the same image, with more analysis points attributed to the wider component.Components with greater intensity also project more points in Fourier space, which can generate anisotropy in Fourier space depending on the difference between intensities (greater intensity differences creating greater anisotropy).Gaps, or other breaks in normally continuous features, play a large part in lowering the order coefficient, with more discontinuities distributing analysis points in a greater number of analysis windows.However, when patterns of "long line" and "short line" features meet, the "long line" pattern has a greater weight in order coefficient analysis.A "wavy" pattern also produces a lower order coefficient in comparison to a perfectly straight pattern.
In terms of stromal collagen imaging, the widths of the collagen fibers will not substantially change between layers or even in the same image.Thus most anisotropy will not be generated from varying fiber width.Fibers on plane will contribute more analysis points due to greater second harmonic signal intensity.Gaps in the fibers introduced during image collection play a sizeable role when calculating the order coefficient, longer continuous fibers weighing more heavily in the calculation of the order coefficient.As the collagen fibers of the cornea relax after extraction from the ocular globe, "waviness" plays a large role in the final value of the order coefficient.

Fig. 2 .
Fig.2.Image Processing Procedure.The edges of the collagen fibers in the raw images (a) were sharpened (b) and a Hanning window (c) was applied to the image to produce a windowed image (d), which concentrates analysis on the central fibers.The windowed image was Fourier transformed in two dimensions (e) and the resulting transformed image was translated to a polar representation (f), where anti-biasing windows were applied (points within 90° ± 2.8° and 270° ± 2.8° on the polar plot were excluded).The polar plot was then divided into 10° windows (g) (30° shown in (g) for easier visibility to the reader) and the points in each window were used to calculate the order coefficient (h).The white bars in (a), (b), and (d) represent a distance of 60 μm.

Fig. 5 .
Fig. 5. Images of collagen from the anterior ((a) and (b)), intermediate (c), and posterior (d) of a de-epithelialized eye.The anterior images were collected at 10 μm and 140 μm below that anterior stromal surface, the intermediate 740 μm, and the posterior 1015 μm.White bars represent 50 μm.All images presented here are from an anterior to posterior z-scan and the contrast was artificially increased for the images presented here to increase fiber visibility for the reader.

Fig. 6 .
Fig. 6.SHG images of stromal collagen fibers ((a), (c), and (e)) with different treatments and corresponding Fourier transforms ((b), (d), and (f)).The untreated ((a) and (b)) and riboflavin only ((c) and (d)) treatments showed more fibers with crimps and the Fourier transform was spread across more windows than the UVX case ((e) and (f)).The order coefficient values of the individual images were; 0.357 for the untreated, 0.369 for the riboflavin only, and 0.791 for the UVX cornea.Images have been rotated by 90° and contrast has been enhanced so the reader can more clearly see the collagen fibers.All three images were taken 108 μm below the anterior surface of the stroma.The white bars represent a distance of 50 μm.

Fig. 9 .
Fig. 9. Order coefficient averages of corneal treatments.The number of images used to compute each average in the separate depths regions and treatments (averages computed across all samples in each treatment, untreated n = 8, riboflavin only n = 8, UVX n = 8) is listed in Table 1.Black bars represent the variance.

Fig. 10 .
Fig. 10.Ideal examples of order and disorder with Fourier transforms.

Fig. 10 (
Fig.10(o) is an image of the one-pixel width black and white bars with short (2-5 pixel) intensity blocks overlaid on a sinusoidal pattern with amplitude of 5 and wavelength of 60 pixels.The corresponding Fourier transform is shown in Fig.10(r) and the corresponding order coefficient is 0.378."Waviness", as demonstrated in this example, plays a noticeable part in lowering the value of the order coefficient.To review, the order coefficient has the lowest value possible when the distribution of analysis points is evenly distributed and rises when anisotropies are present.The order