Analysis of spatial lamellar distribution from adaptive-optics second harmonic generation corneal images.

The spatial organization of stromal collagen of ex-vivo corneas has been quantified in adaptive-optics second harmonic generation (SHG) images by means of an optimized Fourier transform (FT) based analysis. At a particular depth location, adjacent lamellae often present similar orientations and run parallel to the corneal surface. However this pattern might be combined with interweaved collagen bundles leading to crosshatched structures with different orientations. The procedure here reported provides us with both principal and crosshatched angles. This is also able to automatically distinguish a random distribution from a cross-shaped one, since it uses the ratio of the axes lengths of the best-fitted ellipse of the FT data as an auxiliary parameter. The technique has successfully been applied to SHG images of healthy corneas (both stroma and Bowman's layer) of different species and to corneas undergoing cross-linking treatment.

1. Introduction

Collagen (type I) is the main structural element of the corneal stroma. It presents a natural non-centrosymmetric organization with a triple helix fibril structure which generates a strong second harmonic generation (SHG) signal [1]. These fibrils assemble to form collagen bundles or lamellae [2]. SHG microscopy enables imaging of corneal collagen structures without staining procedures [1,3-6]. Corneal lamellae are often organized and their orientations are often analyzed in a qualitative manner. However, quantification facilitates extraction of sensitive changes in collagen fiber organization due to pathology [4,7], surgery [8-11] or damage [12,13], and, thus may help in early diagnosis and follow-up processes.

There has been an increasing interest in quantifying corneal lamellae organization from SHG imaging under different experimental conditions, due to possible clinical applications. Thermal responses in corneal SHG images have been characterized by means of a fast Fourier transform (FT) algorithm [12]. The aspect ratio (AR) computed as the quotient between the short and long axes of the ellipse fitting the FT spectrum image, was used as a parameter to compute changes in lamellar organization. As expected, the features in the spectral images moved from an ellipse to a circle as the collagen distribution tended to present a random pattern. Rao and associates used a simpler method: the best-fit line to the binary spectral image (which is perpendicular to the collagen preferred orientation) [14]. As indicator of the number of fibers that deviate from the preferred orientation they used the standard fitting error. A FT analysis was also recently made to compare structural alteration between normal and keratoconic corneas imaged with SHG microscopy [15]. The AR served as a quantitative measure of fiber direction determination.

Most of these methods combined analytical and manual algorithms, which might be misleading or bias results if different operators are involved. Moreover, the use of the AR parameter often fails when the collagen fibers are highly crimped (wavy). In that sense, a changeling method to overpass this limitation has recently been reported [16]. This technique is based on the Radon transform of the FT image.

Recent experiments based on SHG imaging have reported that healthy corneal collagen is not always arranged in a well-organized structure, but in more complex and heterogeneous patterns, including interweaving, branching, undulations (“crimps”) and crosshatching, among others [6,15,17,18]. These features are better visualized when combining adaptive optics (AO) or wavefront control techniques with multiphoton microscopy [6,19]. Moreover, the use of the AR as previously reported [12,15] might erroneously identify the structure as a random distribution in corneas presenting crosshatched structures.

In this sense, the present work goes a step further into the analysis and quantification of collagen fiber organization in SHG imaged corneal tissues. We propose an optimized and refined method based on a FT-based algorithm. This procedure will enable computation of the preferential direction of corneal lamellae, as well as explore the presence of crosshatched structures and distinguish them from random ones.

2. Experimental procedure

A research AO multiphoton microscope was used to acquire SHG images of ex-vivo corneas of different species (humans included) in the backward direction [19]. The system was completely automated and controlled through custom software. The laser average power for SHG imaging depended on the sample and ranged between 100 and 130 mW. The use of AO enables imaging of deeper layers (posterior stroma) with enough contrast to be processed. An example showing the benefits of using AO in SHG corneal microscopy is presented in Fig. 1.

Fig. 1

Effects of using AO correction to improve the visibility of features within the corneal stroma (eagle). AO off (a); AO on (b). Images correspond to a depth location of 200 μm.

Details on the experimental system, procedure and tissue maintenance can be found elsewhere [6,11,19]. In brief, porcine and bovine corneas were obtained from a local slaughterhouse. Intact eyes were placed upside down on a glass bottom dish and filled with a solution combining Hanks salts, sodium and L-glutamine. Corneas from rats, rabbits, and chickens were excised with a trephine and immersed in the same solution. Normal human corneas from donors not suitable for transplantation were kindly provided by the “Servicio de Oftalmología,” Hospital Universitario Virgen de la Arrixaca, Murcia, Spain. After trephination, samples were stored in Optisol solution.

All imaged areas within this work are 180x180 μm in size, correspond to the central part of the cornea (although different depth positions) and were acquired with the AO module in operation (Fig. 1(b)). Each initial SHG image involved in the procedure described below was the average of three individual frames. The entire image processing was performed with MatLabTM (The MathWorks, Inc, Natick, MA). The optimized procedure developed to analyze the distribution of the corneal stroma is described in the following.

The initial averaged image (see Fig. 2(a)) was first filtered by a Gaussian low-pass filter to reduce background noise. A Hamming window was applied to the image prior to Fourier analysis in order to avoid artifacts. This image (the entire one or, alternatively, any specific region of interest or equally sized domains) was then subjected to two-dimensional FT. For each image the logarithm of the absolute value of the FT image was calculated and visualized (Fig. 2(b)).We then proceeded by thresholding the FT image to a binary matrix, centered at (0,0) and plotted in polar coordinates (r, θ) (Fig. 2(c)).

Fig. 2

Example of the procedure steps. (a) SHG image (bovine cornea, 330-μm depth); (b) FT spectrum image; (c) thresholded and binarized FT image (blue dots) with the best fitted ellipse (dashed black line), the scanned window (red lines) and the computed axes (green lines).

Due to the symmetry of the FT, for the range θ = [0°, 210°], the binary data points were scanned with a 30°-range window (Fig. 2(c)), moving at step sizes of 1°, counting the number of points in the window every step of the way. If a preferential direction exists in the FFT image, the plot of the number of points for each scanned window as a function of polar angle presents a maximum for that direction. Once the window with the highest number of points was found, the angle of the best-fit is the average of the angular location of ten points with highest r (i.e. highest frequency points). This number of points was determined after some preliminary tests with a set of corneal tissues under different experimental conditions providing different collagen patterns (healthy, edematous, dehydrated). This averaged angle θFT between the positive X-axis and the line of best-fit lies 90° apart from the preferential direction (angle φpr) of the corneal lamellae within the area under analysis. Once φpr is determined, the AR of the best fitted ellipse is computed. Taking previous published data into consideration [12,15], if the computed AR value was smaller than 0.50, φpr was chosen as the preferential direction of the lamellae and the analysis of the SHG image finishes.

Alternatively, if AR was higher than 0.5, we proceed to find the line of best fit for crosshatched structures (if any, see for instance Figs. 2(a), 3(b) and 4(a)). Since crosshatched structures might not be exactly at 90° from φpr, the approach was similar. The only change is the range of data over which the counting window scans. A forbidden window of ± 30° is then defined around the best-fit angle computed before. A 30°-range window will again scan and count points across the rest of the 180°-field. If within the final chosen window all the points with highest frequencies are outside the best fitted ellipse, the averaged angle θFT is computed again. The angle of crosshatching (φcrh) lies 90° apart from this angle. If a crosshatched structure is present, the number of points for each scanned window must have a secondary maximum around that direction. Otherwise the distribution of corneal collagen is considered as random and the computed value for φpr is not taken into account. For this case all frequencies in the FT image are equally distributed (see Fig. 3(f) as an example).

Fig. 3

SHG images of corneal collagen in porcine (a), human (b) and rat (c) eyes (depth locations: 150, 275 and 100 μm). Bottom panels are the corresponding FT images (see text).

Fig. 4

SHG microscopy image of (a) corneal stroma (chicken, 70 μm-depth); (b) Bowman’s layer (rabbit); (c) stroma after CXL (porcine, 150-μm depth). Bottom panels are the corresponding 2D FT maps computed using the procedure here reported.

3. Results

In order to gain further insight into the stromal arrangement, a quantitative analysis of lamellar distribution for corneal tissues of different species and experimental conditions was carried out. The steps of the proposed method (Fig. 2) were then used to compute the dominant angle of stromal lamellae, as well as the direction of crosshatching (if any).

Figure 2(a) presents a SHG corneal image from a bovine. Collagen fibers running close to both the horizontal and vertical directions can be observed. This simple observation indicates that both principal and crosshatched axes might be present. These two directions are clearly present in the FT image (Figs. 2(b) and 2(c)). A value of φpr = 102° was obtained following the procedure described above. Since AR parameter was 0.72 (see ellipse in Fig. 2(c)), the customized software proceeded to achieve the angular information for the crosshatched structures. For this particular case, φcrh~178° (or −2°). Axes (green lines) depicted in Fig. 2(c) are 90° apart from φpr and φcrh.

Figure 3 shows SHG images of corneal stroma in three different species, together with the corresponding 2-dimensional FT images. The regular FT analysis based on computing the AR of the best fitted ellipse provides a value of 0.15 for Fig. 3(a) and a dominant direction of the lamellae of ~30° (90° apart from the set of points of the FT image of Fig. 3(d)). Since AR = 0.15, the software automatically stops and provides a unique angle φpr.

On the contrary, AR for the SHG image in Fig. 3(b) gives a value of 0.79 (i.e. an ellipse closer to a circle), what would mean that the distribution of collagen for this cornea is rather random [12]. However it can be observed that the corresponding FT image (Fig. 3(e)) shows two preferential directions (crosshatched pattern) at 111° and ~33° due to the crosshatched appearance of collagen fibers. This is a misleading interpretation originating from having a fitted ellipse in the FT image with long and short axes of similar lengths.

For the sense of completeness Fig. 3(c) presents a cornea in which collagen distribution does not show any dominant direction. As expected, the scanning window used in the procedure did not find a location with a maximum number of points (Fig. 3(f)). This was confirmed by means of the AR value directly computed from the FT image (0.95).

Figure 4(a) shows another example of stromal distribution containing interweaving (a chicken cornea). Automatically computed values for AR, φpr and φcrh were 0.73, ~142° and ~44° respectively. Figure 4(b) presents a SHG image corresponding to the Bowman’s layer of a rabbit cornea. This layer is located between the corneal epithelium and the anterior stroma and is made of short fibrils with random orientation and a great degree of interweaving [2,4,6]. A simple visual inspection shows the absence of a preferential collagen direction. This is clearly shown in the shape of the FT map where AR = 0.79 (neither φpr nor φcrh were computed).

Finally, Fig. 4(c) depicts the collagen distributions for a cornea after cross-linking treatment (CXL) [10,11]. CXL procedure is successfully used to stabilize keratoconus progression by increasing corneal stiffness [20]. Noticeable changes appear in the arrangement of the collagen fibers after CXL procedure when compared with the regular orientation of an untreated cornea in Fig. 3(a). Although a predominant direction still exists (for this case horizontal in the SHG image and vertical in the FT one), the map of fibers presents lower spatial frequency (i.e. decrease in the packing of the fibrils). This can be seen in Fig. 4(f). For this particular specimen AR = 0.45 and φcrh~172° (or −8°). Since AR<0.5, the computation of a crosshatched direction is omitted.

4. Discussion and conclusions

An optimized algorithm has been developed to quantitatively analyze the stromal organization of SHG cornea images acquired with an AO multiphoton microscope. The procedure is based on a FT method. This includes not only the AR of the best-fit ellipse, but also the computation of collagen preferential directions in crosshatched structures by means of a scanning window covering the FT image. Unlike previous algorithms, the present one is able to automatically distinguish between random distributed structures and those presenting interwoven or crosshatched distributions.

The AR has often been used as a parameter to quantify corneal collagen architecture. However this work has showed that AR itself might lead to erroneous interpretations (Figs. 2(a), 3(b) and 4(a)). Apart from AR, other quantitative markers (also computed from FT images) such as entropy or disorganization parameter [12,21] and the maximum spatial frequency [14] have been reported to quantify lamellar organization.

Although the procedure has been applied to images of the same size (180x180 μm), our custom software allows the operator to create multiple windows with different sizes and move them across the image in order to run the procedure at any region of interest within the image, or cover the entire image with multiple domains. Smaller windows allow the determination of local fiber directions and are appropriate for samples containing neighbouring areas with different orientations. Moreover, since the anterior stroma and the Bowman’s membrane present very different collagen distributions this method can be used to accurately determine the thickness of this corneal structure by determining a border between ordered and random layers.

The interest in collagen distribution has not been restricted to the cornea. Other collagen-based tissues with biological or clinical interest were also investigated. These include dermis connective tissue [22,23], tendons [24-26], ear cartilage [24], sclera [24,27] and breast tissue [28] among others. Most of these studies used the FT to extract numerical information from collagen fiber organization. Alternative parameters such as the orientation ratio index [22], the FT image covariance matrix [23], the power spectral gradient [25] or the degree of FT eccentricity [27] were also used. Although the present algorithm has been applied to corneal tissues, its implementation in the study of these collagen-based tissues might also be possible. Since these tissues could combine both well-ordered structures and random distributions, the use of the reported algorithm would help to improve the quantitative analysis. However, this FT procedure is not only restricted to SHG imaging. It can also be used with any type of biological tissue or material imaging if the calculation of preferential orientations is required. Moreover, applications to distinguish healthy from pathological tissues might also be of interest.

It has been reported that despite corneal transparency, the lamellar distribution might not be as regular as expected [6]. Moreover the alignment of lamellae is particular for different species [6] and can be modified by different factors including surgery [11], pathologies [4,15], thermal damage [12,13], edema [29], intraocular pressure changes [30] and biomechanical conditions [18]. In this sense, an optimized quantitative method as the one here reported could help to investigate differences in stromal collagen fibers under a variety of experimental conditions.

Some biological tissues show wavy structures and these might also appear in pathological specimens. The reported procedure could also deal with the analysis of these particular patterns. If a wavy arrangement is present, the analysis of the image must be “zonal”, since a regular “global” analysis might lead to erroneous conclusions. With this local analysis multiple windows across the image are created to define adjacent sub-images. This will allow a more appropriate and accurate study of images containing waving structures.

In conclusion, the procedure described here provides a useful and reliable tool for characterizing corneal stroma organization. A special emphasis was put in differentiating random distributions from crosshatched arrangements (and determining the collagen preferential directions). Lamellar distribution was effectively quantified for both normal and surgically-altered specimens. Our results suggest that the parameter AR directly derived from FT analysis of SHG corneal images might fail as an indicator of stromal distribution. AR is also computed in our procedure; however, this is only used as an indicator of a possible existence of interweaving structures. SHG images here shown mostly correspond to intermediate corneal locations (central stroma). AO allows the acquisition of SHG images from deeper layers with enough contrast to be processed [31]. Then this algorithm could be used with plane-by-plane SHG images covering the entire cornea. This would provide helpful information on the collagen distribution as a function of depth. The use of this method in clinical environments might also be of special interest for quantifying changes in lamellar organization after surgery (as shown for CXL), accidental injury or damage, and disease (keratoconus, edema).