1. Introduction

As of 2010, there were an estimated 285 million visually impaired and 39 million blind individuals worldwide [1]. In the United States, proliferative diabetic retinopathy (PDR) and wet age-related macular degeneration (AMD) are two of the leading causes of severe vision-loss and blindness [1, 2]. PDR is a complication of diabetes that is characterized by neovascularization originating from the retina and optic disc, which can result in hemorrhage, fibrosis, traction, and retinal detachment [3, 4]. An estimated 415 million adults suffer from diabetes globally, and almost half of all diabetics are expected to experience some degree of retinopathy [4, 5]. AMD is the leading cause of blindness among people aged 55 and older in the developed world and affect more than 1.75 million individuals in the United States [6, 7]. Wet AMD is characterized by neovascularization, which often leads to hemorrhage and exudation. While only 10 percent of AMD patients develop wet AMD, severe vision-loss progresses quickly in the majority of these patients [6].

Vascular endothelial growth factor (VEGF) inhibitors have become standard treatments for both PDR and wet AMD [8, 9]. Intravitreal anti-VEGF injection has been shown to maintain vision, with 91.5-95.4% of wet AMD patients in one clinical study showing less than a 15-letter decrease in visual acuity (VA) from baseline [10]. However, long-term anti-VEGF treatment has been suggested to be less effective, and studies have shown that more than 30% of patients experience a VA loss of >15-letters after more than five years of treatment [11, 12]. Another major limitation of anti-VEGF therapy is the need for repeated injections on either a monthly or as-needed basis (~6-7 injections per year) [4, 10, 13]. In addition to patient anxiety, repeat injections also increase the risk of side effects such as endophthalmitis [14], acute increases in intraocular pressure requiring topical or surgical anti-glaucoma interventions [15], and off-target drug effects including loss of retinal ganglion cells and circulation disturbances in the choriocapillaris [9]. The aforementioned limitations of the current clinical standard-of-care and lack of understanding of the structural, metabolic, and vascular changes underlying retinal neovascularization highlights the need to identify mechanisms of pathogenesis and novel anti-angiogenic therapies.

The zebrafish (Danio rerio) is a popular model organism because mutant phenotypes of many human diseases are readily available and its fecundity and life cycle make it well-suited for large-scale studies [16, 17]. As an ophthalmological model, the zebrafish retina shares similar structure and function with that of humans and other vertebrates [18]. Similar to humans, the zebrafish retina is composed of seven major cell types (Müller glia and six neural cell types), three nuclear layers separated by two plexiform layers, and a highly ordered mosaic of neurons in each layer [19, 20]. As diurnal species, humans and zebrafish both have cone-dominant vision, in contrast to rod-dominant vision found in mice, which has advantages for studying cone degeneration diseases such as AMD [20, 21].

Zebrafish are also affected by and readily absorb compounds from their aqueous environment, which enables induction of pathologies and delivery of chemical compounds without the need for injections [22]. Retinal vasculopathies can be modeled by exposing animals to hypoxic water for 3-10 days to induce neovascularization and vascular leakage [21]. Similarly, exposure to glucose induces hyperglycemia, which has been shown to result in retinal structure abnormalities similar to those in DR [21]. As a pharmacological model, 82% of disease-causing human proteins have orthologues in zebrafish, and pharmacological effects are highly conserved between humans and zebrafish [23–25].

The majority of drug screens in zebrafish are performed using larval animals because their transparency and size are well-suited for imaging and housing in large-scale studies [26]. However, normal vascular development [27] and differences in inflammatory and immune responses between larva and adults [27] may confound structural and functional changes in larval zebrafish models of neovascularization. In this study, we demonstrate in vivo retinal imaging in adult zebrafish (≥3 months post-fertilization, mpf) using optical coherence tomography (OCT) [28] and OCT angiography (OCT-A) [29], and present post-processing algorithms for vascular segmentation and biometry. Quantitative measurements of retinal perfusion and angiogenesis during longitudinal studies can provide insights into disease pathogenesis and therapeutic efficacy in drug screens for novel anti-angiogenic compounds. In addition to tracking functional changes, retinal vascular biometry can also be used as a method for uniquely identifying individual animals with high sensitivity and specificity. To this end, we developed and validated a novel identification method that obviates the need for physical marking methods such as elastomer marking, freeze branding, removal of specific scales, fin clipping, and dorsal fin tagging [30]. We believe the retinal vascular biometry methods presented here are robust enabling technologies that will broadly benefit large-scale zebrafish studies.

2. Methods

2.1 Imaging system

OCT and OCT-A imaging of wild-type zebrafish was performed using a custom-built spectral domain OCT (SD-OCT) system (Fig. 1(a)). A superluminescent diode (InPhenix) with 855 nm central wavelength and 90 nm bandwidth at full width half maximum was split between reference and sample arms using an 80:20 coupler, respectively, and detected on the central 2048 pix. of a 4096 pix. linear CMOS array with 125 kHz line-rate (spL4096-140km, Basler). Measured system SNR was 91.7 dB with 700 µW at the sample. Measured axial resolution was 2.56 µm (tissue) and −6 dB SNR falloff was at 1.1 mm. Scanning relay optics in the sample arm were designed for a 1 mm diameter spot at the pupil. This includes a 4x demagnifying relay using a 1-inch diameter achromat (f_s) and two 1/2-inch diameter achromats (f_oph) (Fig. 1(a)). Lateral resolution was estimated as 2.8 µm using a model eye but limited by aberrations.

 

Fig. 1 In vivo retinal OCT imaging in zebrafish. (a) Custom-built SD-OCT system. CMOS, detector; f, collimating, objective, ophthalmic, and scan lenses; G, galvanometers; M, mirror; PC, polarization controller; SLD, superluminescent diode; VPHG, volume phase holograph grating. (b) Zebrafish retina was imaged in air through a contact lens and positioned using a custom holder (inset) on a 5-axis alignment stage.

Download Full Size | PDF

2.2 Imaging protocol

In vivo imaging was performed under an animal protocol approved by the Institutional Animal Care and Use Committee (IACUC) at Cleveland Clinic. Both eyes in ten adult wild-type zebrafish (≥3 mpf) were imaged repeatedly on ten different days over four weeks (10 total data sets per eye). Repeat imaging days were separated by 48 hours. During each imaging session, animals were anesthetized using a 0.14% Tricaine solution. Retinal imaging was performed in air and a polymethyl methacrylate contact lens (Cantor & Nissell, Ltd., Northamptonshire, U.K.) was used to compensate for optical power differences in the zebrafish eye when imaging in air instead of water (Fig. 1(b), inset) [31]. Zebrafish OCT and OCT-A volumetric data sets were centered on the optical nerve head (ONH) using a 5-axis stage (Fig. 1(b), Leica Microsystems). OCT volumes consisted of 2500 B-scans (2048 x 500 pix.) acquired in approximately 10 s and covered 10° field-of-view (FOV). These data sets included five repeated B-scans at each lateral position for OCT-A post-processing. Animal imaging and handling were performed in less than 10 minutes followed by recovery, during which water was forced across the gills to ensure animal survival for the duration of the study [31]. Between imaging days, fish were housed in light and temperature controlled rooms in separated boxes and grouped so that individual animals were uniquely identifiable by strip/spot patterns and caudal fin cuts.

2.3 Vessel segmentation and labeling

Vessel segmentation, labeling, and feature extraction were performed on en face OCT-A projections using custom developed algorithms (Fig. 2). The steps are grouped into four categories to better represent the processing pipeline. Each of these steps are described in detail in the following sections.

 

Fig. 2 Block diagram showing vessel segmentation, labeling, and feature extraction steps for automated vascular biometry.

Download Full Size | PDF

Data set acquisition and processing

Bulk motion in OCT B-scan were removed using a discrete Fourier transform cross-correlation [32]. First, repeated B-scans at the same location were registered to each other and averaged. Then, all averaged B-scans in the volumetric data set were registered to adjacent scans. Volume registration parameters were calculated using OCT data and applied to the corresponding OCT-A frames.

OCT-A vasculature maps were calculated using weighted optical microangiography (wOMAG) [33]. wOMAG uses OCT intensities to remove OCT-A artifacts caused by tissue hyper-reflectivity. Here, optical microangiography frames were weighted by an intensity decorrelation function, (D/D₀)ⁿ. The decorrelation coefficient, D, represents the intensity decorrelation between repeated B-scans, and optimal values for D₀ and n were found experimentally (D₀ = 0.1, n = 0.5).

Skeletonization of vascular maps

En face OCT-A projections were lowpass filtered to smooth vessel contours and remove speckle noise. A vertical intensity gradient was also calculated and subtracted from the filtered vessel maps to remove breathing artifacts (horizontal streaks). The resulting OCT projections were then binarized and skeletonized using dilation and erosion [34]. Morphological dilation and erosion expand and compress binary images. OCT-A projections were dilated using a circular kernel (5 pixel radius) to remove vessel thickness nonuniformity and then eroded to obtain a one pixel width skeleton.

Vessel classification and branch point detection

Skeletonized vasculature maps were used to detect vessel branch points. Vessels were first classified as either ONH or retinal vessels. The retinal pigment epithelium (RPE) was segmented from OCT cross-sections [35] and the ONH was identified by discontinuities in the RPE. The resulting ONH segmentation was then fit to a circle that defined the ONH center and radius. All vessel branches within the ONH radius were classified as ONH vessels and all remaining vessels were further processed for branch point identification.

A set of 18 predefined 3x3 pixel branch kernels was created to represent all possible orientations of vessel bifurcation and trifurcation (Fig. 3). The spatial location of each vessel branch point was identified by convolving the skeletonized vasculature map with each kernel. Branch points were then used to classify each vessel segment by branch generation relative to the ONH (Fig. 4).

 

Fig. 3 Branch kernel matrix. 18 predefined 3x3 pixel kernels were created to represent all possible vessel branch orientations.

Download Full Size | PDF

 

Fig. 4 Quantitative vascular biometry. Length, curvature, and angle were defined using skeletonized vascular maps. Branch generations were assigned based on their connectivity relative to the ONH.

Download Full Size | PDF

Quantitative vascular biometry

Quantitative biometry was performed by extracting vessel segment length, curvature, and angle between branch points in the skeletonized vessel map (Fig. 4). Segment length was defined as the total number of pixels in each vessel segment, curvature was calculated as the ratio between the vessel segment length and Euclidean distance between corresponding branch points, and angle was calculated as the angle between vessel segments. Vascular biometrics were extracted for each vessel branch originating from the ONH starting from the 12 o’clock position and moving clockwise around the ONH.

2.4 Retinal vascular biometry for animal identification

Retinal vascular biometry was evaluated as a robust, noncontact, and noninvasive method for unique zebrafish identification over longitudinal timepoints. A Pearson correlation coefficient matrix was calculated by comparing vessel branch length, curvature, and angle for all branch generations between all data sets. We then implemented different methods for comparing vascular biometry data between branch generations to uniquely identify eyes across our longitudinal data sets.

A weighted average method was implemented to combine correlation values for each generation into a single correlation matrix between all data sets. Generation weights were calculated by optimizing identification sensitivity and specificity using particle swarm optimization [36] and a genetic algorithm [37]. Both methods found 0.9852, 0.0214, and 0.0001 as the optimal weights for first three generations, respectively, and zero for all higher generation branches. Here, the weighted average favored contributions from lower generations to compensate for variability in OCT/OCT-A FOV between longitudinal timepoints, which may result in inconsistent biometrics from higher generation vessel branches.

In addition to the weighted average approach, we also explored other methods for comparing correlation values that did not require a training set for branch generation weight optimization. The biometric identification accuracy of the weighted average method depended entirely on the robustness of the training set and can be easily biased. Based on our weight optimization results, it was clear that biometric identification should be dominated by first generation branches. Similarly, when comparing longitudinal data sets, the highest correlation values from matching eyes were also found in the first-generation branches. Thus, we simply used the maximum correlation coefficient from the first three branch generations as the correlation coefficient for the corresponding data set.

For both the weighted average and maximum correlation approaches, a combined correlation matrix comparing all 200 data sets was created (i.e., 200 x 200 matrix). Using a priori knowledge that imaging was repeated in each eye ten times, the ten highest correlation matrix values for each data set were selected to correspond to matching eyes. Sensitivity, specificity, and accuracy were then calculated to assess the performance of our biometric identification algorithms.

3. Results

A representative in vivo OCT data set centered on the ONH is shown in Fig. 5 with retinal layers labeled based on previous studies comparing OCT cross-sections to corresponding histology (Fig. 5(b)) [38, 39].

 

Fig. 5 In vivo retinal OCT in zebrafish. (a) En face OCT projection with representative orthogonal cross-sections (blue/red lines and insets). (b) 5-frame averaged OCT cross-section with labeled retinal layers. GCL, ganglion cell layer; IPL, inner plexiform layer; INL, inner nuclear layer; OPL, outer plexiform layer; ONL, outer nuclear layer; OS, outer segment; and RPE, retinal pigment epithelium.

Download Full Size | PDF

The corresponding en face OCT-A projection shows central major vessels radiating outward from the ONH (Fig. 6). OCT-A B-scans show vessel cross-sections at the surface of the retina and flow artifacts at the RPE below each vessel. We distinguish retinal vessels from artifacts by segmenting and isolating only the upper layer of flow signals (Fig. 6(a), arrows) because in the adult zebrafish, the retinal vasculature forms a membranous layer that is attached to the vitreal interface [40].

 

Fig. 6 In vivo retinal OCT-A in zebrafish. (a) En face OCT-A projection with representative orthogonal cross-sections (blue/red lines and insets) showing retinal vessels (arrow) and RPE artifacts. (b) OCT-A projection with corresponding segmentation mask overlay. Vessel branches are color-coded based on their branch generation relative to the ONH (white circle). Scale bar 1°.

Download Full Size | PDF

OCT-A vessel maps were skeletonized, and the resulting vessel segments were color-coded to represent different branch generations relative to the ONH (Fig. 6(b)). Any skeletonized vessel segments that began and ended within the ONH (white circle) were ignored. Vessels that began inside and ended outside of the ONH were considered the first branch generation (green) and branch generations were incremented radially outward from the ONH. Total processing time for each OCT-A data set was approximately 13 minutes, the bulk of which was spent on volumetric registration (>12 min.). Skeletonization and vessel segmentation was performed in ~30 s per vascular map.

Automatically segmented vessel maps were evaluated by manual graders to quantify the robustness of our algorithm. Errors were identified in 2.5% of data (5 of 200 vessel maps) and classified as either segmentation errors at the periphery of the FOV (Fig. 7) or inside the ONH (Fig. 8). At the edge of the FOV, segmentation is confounded by areas of low OCT-A contrast and cropped vessel branches. Poor contrast results in missed branches and branch points (Fig. 7, pink arrows). Similarly, vessels cropped at the end of the imaging FOV may result in missed branch points (Fig. 7, orange arrows). Missing branches or branch points confound all downstream analyses, including branch generation labeling and quantification of branch length, curvature, and angle. Mislabeled branch generations also occurred because of overlap between a vessel branch point and the ONH rim. In these cases, there was ambiguity in classifying the vessel segment as either a retinal or ONH vessel, which led to the mislabeling of subsequent branch generations originating from the initial mislabeled ONH branch (Fig. 8).

 

Fig. 7 Segmentation errors at the FOV periphery. Representative en face (a), (d) OCT; (b), (e) OCT-A; and (c), (f) segmentation maps of the same eye at two longitudinal timepoints. Some generation 3 branches were misidentified as part of generation 2 branches because of poor OCT-A contrast (pink arrows). Similarly, branches cropped by the edge of the FOV (dashed lines) with sufficient contrast may be misidentified as part of the preceding generation (orange arrows).

Download Full Size | PDF

 

Fig. 8 Segmentation errors at the ONH. (a)-(c) En face OCT-A, (d)-(f) segmentation maps, and (g)-(i) magnified views around the ONH of the same eye at three longitudinal timepoints. (d), (g) Correctly labeled data sets show ONH vessels in gray. (e), (h) and (f), (i) Branch generation labeling errors show ONH vessels labeled as first-generation vessels (arrows).

Download Full Size | PDF

Automatically extracted quantitative vascular biometrics were color-coded and plotted as maps of vessel branch length, curvature, and angle to enable qualitative comparisons between longitudinal data across different animals and eyes (Fig. 9). Biometrics were grouped by left (OS) or right (OD) eye in each fish, and all 10 repeated longitudinal data sets in each eye were displayed as vertical columns. In the horizontal axis, branch segments (row) were grouped by branch generations. Within each eye, similar biometric patterns were observed in repeat data sets, but data between eyes and fish showed significantly different biometric features.

 

Fig. 9 Quantitative vascular biometry. Vessel branch (a) length, (b) curvature, and (c) angle of 10 repeated longitudinal data sets (columns) for both eyes of 10 zebrafish (column groups). Each column group showed biometric changes in 10 imaging sessions for one eye. Biometrics in each eye were grouped by branch generation (row groups) and each row in generation represent one vessel in the respective generation. Color coded biometric features shows change in length, curvature, and angle in each vessel over time.

Download Full Size | PDF

Branch vessel length was identified as the most robust biometric parameter for comparing vascular data between longitudinal timepoints because it had the highest dynamic range and lowest noise of all three parameters (Fig. 9.). Qualitatively and expectedly, biometric data becomes noisier in higher generation (>3) vessel branches because of aforementioned segmentation errors at the edge of the FOV (Fig. 10).

 

Fig. 10 Vessel branch length comparison in the same eye between longitudinal timepoints. (a) Representative biometric map showing vessel branch length differences (columns) between 10 repeated data sets (rows) in one eye. (b)-(d) Segmented vascular maps at 3 representative timepoints showing similar length parameters in the first and second vessel generations and increased variability in higher generations as a result of FOV differences (arrow and asterisks).

Download Full Size | PDF

Higher branch generation biometric parameters were strongly influenced by eye position in the FOV. As explained in Section 2.4, the first-generation branches dominated biometric identification as compared to higher generations. Our weighted average approach achieved 96.5% sensitivity and 99.8% specificity for uniquely identifying data sets from the same eye across longitudinal timepoints. When the previously described segmentation errors were corrected, the weighted average sensitivity and specificity increased to 98.4% and 99.9%, respectively. Our maximum correlation approach achieved 98.9% and 99.9% sensitivity and specificity, respectively, and 99.9% accuracy using the same set of data without the need for a training set. More in-depth analysis of the maximum correlation approach showed that first generation branch correlations were dominant in 87% of all correctly identified data sets and the second and third generation correlations were dominant in 12.6%, and 0.4%, respectively. Similarly, first generation branch correlations were dominant in all misidentified data sets.

4. Discussion

We presented a novel zebrafish identification method that is efficient and robust. Total processing time for vascular segmentation and labeling was ~30 s, and our maximum correlation identification algorithm achieved a sensitivity and specificity of 98.9% and 99.9%, respectively, without the need for a training set. Our proposed method addresses major limitations in large population imaging studies of adult zebrafish, specifically the need to uniquely identify animals between longitudinal timepoints. High sensitivity and specificity identification methods may eliminate the need for fin cutting and housing animals in limited groups (4-6 fish). OCT-A also provided functional information to enable quantitative biometrics in zebrafish models of retinal vascular pathologies. The full retinal vasculature is visualized on OCT-A. Zebrafish retinal vasculature forms a single layer attached to the vitreal interface [41]. Main branches arborize from the ONH, but do not form subretinal plexi as compared to other model animals [31].

In our preliminary study, identification errors were primarily a result of differences in FOV between longitudinal timepoints. Translational differences in FOV affected the visibility of peripheral vessels, which contributed to increased noise in biometric parameters for higher branch generations. Vessel branches that were cropped at the edge of the FOV may cause branch points to be missed during automatic feature identification. The current algorithm defines branch points as the point where a parent branch ramifies into two or three children branches. Loss of children branches by FOV cropping can result in a parent and child branch to be considered as a single vessel branch, resulting in significant errors in biometric quantification (Fig. 10). We overcame these errors by omitting higher generation branch vessels that intersect the FOV edge from our identification algorithm to avoid the potential for associated errors. Rotational differences of the FOV at different timepoints also limited the accuracy of biometric comparisons because branch vessels may be mislabeled between data sets. Rotational errors may be addressed using volumetric registration of images at different timepoints at the expense of increased computational complexity. ONH segmentation errors were also observed between timepoints. Inaccurate identification of ONH vessels resulted in vessel generation labeling errors for all subsequent branch vessels. One possible solution for ONH vessel segmentation errors may be to appropriately threshold first generation branch length because ONH vessels are significantly shorter than first generation vessels.

The major computational bottleneck of our current algorithm was volumetric registration, which was ~12 min per volume. However, registration does not significantly limit the utility of our biometry method because it is performed in post-processing and real-time operation is not required. In addition, we observed that registration of repeated OCT-A frames at the same position was not necessary in all eyes. Thus, increased OCT/OCT-A imaging speeds may obviate the need for volumetric registration to remove bulk motion noise.

5. Conclusion

We demonstrated the utility of OCT and OCT-A for structural and functional imaging of the zebrafish retina. Quantitative vascular biometric parameters were automatically extracted from OCT-A vessel maps and used to uniquely identify eyes in different fish between longitudinal timepoints. Using a maximum correlation approach, we achieved an identification accuracy of 99.9% in a preliminary study of wild-type fish over 10 repeated timepoints. This technology eliminates the need to manually mark and identify animals and provides quantitative metrics for studying functional changes in zebrafish models of retinal pathologies.