Mapping the 3D Connectivity of the Rat Inner Retinal Vascular Network Using OCT Angiography

Abstract

Purpose: The purpose of this study is to demonstrate three-dimensional (3D) graphing based on optical coherence tomography (OCT) angiography for characterization of the inner retinal vascular architecture and determination of its topologic principles.

Methods: Rat eyes (N = 3) were imaged with a 1300-nm spectral/Fourier domain OCT microscope. A topologic model of the inner retinal vascular network was obtained from OCT angiography data using a combination of automated and manually-guided image processing techniques. Using a resistive network model, with experimentally-quantified flow in major retinal vessels near the optic nerve head as boundary conditions, theoretical changes in the distribution of flow induced by vessel dilations were inferred.

Results: A topologically-representative 3D vectorized graph of the inner retinal vasculature, derived from OCT angiography data, is presented. The laminar and compartmental connectivity of the vasculature are characterized. In contrast to sparse connectivity between the superficial vitreal vasculature and capillary plexuses of the inner retina, connectivity between the two capillary plexus layers is dense. Simulated dilation of single arterioles is shown to produce both localized and lamina-specific changes in blood flow, while dilation of capillaries in a given retinal vascular layer is shown to lead to increased total flow in that layer.

Conclusions: Our graphing and modeling data suggest that vascular architecture enables both local and lamina-specific control of blood flow in the inner retina. The imaging, graph analysis, and modeling approach presented here will help provide a detailed characterization of vascular changes in a variety of retinal diseases, both in experimental preclinical models and human subjects.

Inner retinal metabolism requires a continuous supply of oxygen and other nutrients, as well as waste removal, by blood flow through a network of vessels.^1,2 Vascular architecture constrains how blood flow can be controlled to meet metabolic needs in the inner retina.³ Moreover, vascular structure in the retina may aid the diagnosis and management of diseases that involve vascular dysfunction, such as diabetic retinopathy,⁴ glaucoma,^5,6 and cardiovascular disease.^7,8 Thus, comprehensive architectural analysis of the inner retinal vasculature will impact both basic research and clinical practice. 

Optical imaging through the transparent ocular media offers the opportunity to quantitatively assess vasculature noninvasively in vivo. The blood vessels of the inner retina exhibit a characteristic three-dimensional (3D) layered structure, with superficial, intermediate, and deep vessel layers, corresponding to the ganglion cell layer/nerve fiber layer, the inner plexiform layer, and the outer plexiform layer, respectively.⁹ A fundamental limitation of 2D photography-based techniques, such as fluorescein angiography,¹⁰ is the lack of depth discrimination,¹¹ due to overlapping signals from the individual retinal layers and the choroid. 

Optical coherence tomography (OCT)¹² enables 3D optical sectioning of the retina with micron-scale resolution,¹³ employing only intrinsic contrast without the requirement for additional contrast agents. The superior axial resolution of OCT distinguishes the layered architecture of the retina, including the capillary plexuses located in the plexiform layers.¹⁴ Optical coherence tomography angiography techniques,^15–18 based on enhancing contrast arising from motion and scattering, were recently developed for selective imaging of only the “functional” vasculature perfused with moving blood cells. Doppler OCT^19,20 has also been shown to perform quantitative assessment of blood flow.²¹ 

Optical coherence tomography–based studies involving 2D segmentation have proven effective for visualization²² and feature detection,²³ however, 3D analysis is required to assess the true vascular topology, more directly related to function. In this work, we derive a 3D representation, or “graph,” of the vasculature in the rat inner retina based on OCT angiography data. Unlike previous segmentation efforts in the retina, the interconnected nature of the multilayered vascular structure is implicitly incorporated into the derived representations and analysis. Using the graph, in conjunction with experimentally-quantified blood flow data (computed using Doppler OCT methods²⁴) and a simple theoretical network model of blood flow, the capability of blood vessel dilation to produce localized and lamina-specific changes in blood flow is examined. 

Methods

Imaging Protocols

OCT Microscope and Animal Preparation.

The rat eye was imaged with a 1300-nm spectral/Fourier domain OCT microscope²⁵ operating at a speed of 47,000 axial scans per second with a transverse resolution of 7.2 μm and an axial resolution of 3.6 μm.²⁶ The use of 1300-nm light, as opposed to 800-nm light, results in less prominent “tail” artifacts beneath large vessels, thus facilitating image segmentation. Male Sprague-Dawley rats (350–430 g) were anesthetized with 40 to 50 mg/kg/h α-chloralose, infused intravenously, while maintained at 37°C with a homeothermic blanket. Pupils were dilated with 1% tropicamide (Mydriacyl; Alcon Laboratories, Fort Worth, TX, USA), and imaging was performed using a contact lens and goniosol on the cornea. Animals were ventilated with 80% air and 20% oxygen, and arterial blood pressure was monitored through a femoral artery cannula. All procedures were approved by the Institutional Animal Care and Use Committee where these experiments were performed, and adhered to the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research. 

Angiography

The aim of OCT angiography is to enhance the contrast of scattering red blood cells (RBCs) in vessels with respect to the surrounding nonvascular tissue. The complex OCT signal can be considered as the superposition of a static component, a dynamic component, and an additive noise component.²⁷ At any given voxel, these three components are assumed to be independent complex random processes. The power in the dynamic scattering component is related to RBC content; therefore, estimation of this power yields a representation of the angiogram. This can be computed through high-pass filtering the complex OCT signal in time, followed by averaging. Our angiography scanning protocol acquired cross-sectional images with 512 axial (z) scans, repeated twice at each of 5120 transverse (y) locations, over a field-of-view of 1.48 × 1.48 mm².²⁶ The total acquisition time was approximately 2 minutes. 

The high-pass filter was achieved through complex subtraction of repeated frames. The interframe time was approximately 11 milliseconds, sufficient for the dynamic scattering component of the complex OCT signal to decorrelate. Thus, the resultant angiogram was largely insensitive to RBC speed and therefore depended mostly on RBC content.²⁸ Bulk axial motion between frames, which leads to a uniform shift in the phase of the complex OCT signal,²⁹ was fixed by computing and correcting for the phase shift between corresponding A-lines in consecutive frames. The volumetric angiography data from the inner retina was averaged and resampled to 1024 × 1024 × 100 voxels. 

Doppler OCT

For the Doppler OCT scanning protocol used in this study, four volumes (256 cross-sectional images with 4096 axial scans each) were acquired sequentially. Each volume acquisition required approximately 25 to 30 seconds. Doppler shifting causes a pure rotation of the OCT signal in the complex plane, and leads to a frequency shift of the power spectrum. Frequency shifts along the fast axis were first estimated and then converted to velocity axial projections. Velocity axial projections were then used to compute volumetric flow. Because the scanning protocol sampled vessels asynchronously with the heartbeat,²⁵ time-averaged volumetric flow could be obtained using the method of en face integration,²⁴ which eliminates the requirement for explicit determination of vessel angle. The time-averaged flow was measured in individual arteries and veins near the optic nerve head. 

Image Processing

The segmentation and vectorization of OCT angiography data present several unique challenges. In general, blood vessels exhibit a wide range of sizes, curvatures, and geometries.³⁰ Multiple scattering tails from large vessels or shadows from media opacities contribute to variations in signal strength across the imaged field. Optical coherence tomography angiography data is also prone to shot noise, as well as speckle noise³¹ compounded with sample motion, which can lead to imperfect suppression of static tissue in the angiogram. Additionally, underperfused vessel segments may not be detected well in the angiogram. The series of image processing steps described below were specifically tailored to mitigate this set of potential confounds. 

Image Enhancement and Binarization

Image enhancement was performed on the angiogram in order to increase the contrast of the vessels within the image before performing segmentation. In this work, the Frangi filter,³² a widely-adopted multiscale vessel enhancement technique, was applied to the image data. Subsequently, an empirically-determined threshold was applied to the enhanced volume, yielding a binary mask representation of the vasculature, which forms the basis for the further computation of vessel topology. Morphologic opening was then applied to mitigate the effect of noise on the mask. Figure 1 illustrates the main image processing steps involved in obtaining the binary segmented vascular mask. Figure 1A shows a maximum intensity projection of the angiogram image data from the rat inner retina. Figure 1B shows a maximum intensity projection of the enhanced data. Figure 1C shows the binary mask computed via thresholding of the enhanced data. 

Figure 1

 

Segmentation of angiogram image data from the rat inner retina. (A) Maximum intensity projection of angiogram image data; (B) maximum intensity projection of enhanced angiogram image data; (C) 3D rendering of the segmented binary mask.

Skeletonization

For blood vessel networks, skeletonization approaches endeavor to compute 1-voxel thick center lines that are equidistant from the vessel boundaries. Typically, the skeletonization process entails systematically removing voxels from the binary mask, until only a single voxel thick structure remains. The medial axis thinning algorithm described by Lee et al.³³ is one such 3D skeletonization method. A MATLAB (MathWorks, Natick, MA, USA) implementation of this method³⁴ was used to obtain a 3D center-line mask that preserved the connectivity of the initial binary mask. 

Skeletonization errors typically take the form of small erroneous branches occurring due to inhomogeneities in the structure of the binary vessel mask, and gaps in segmented branches (e.g., due to low signal contrast) that are manifested as unconnected endpoints in the skeleton. Many of the small erroneous branches can be removed by simply defining a well-considered minimum length threshold. Threshold relaxation³⁵ and tensor voting³⁶ are two strategies that have been advanced for the correction of vessel-gap errors, though few results have been reported for 3D vascular data.³⁷ Manually-guided approaches allow for comprehensively and accurately correcting skeletons, but are labor-intensive.³⁸ We developed a suite of automatic and manually-guided methods to correct the vascular skeleton, as illustrated in Figure 2. Figure 2A shows uncorrected and corrected skeletons, overlaid on the binary mask. Figure 2B illustrates the resolution of a skeleton error performed by combining multiple skeletonizations at different thresholds. Figure 2C illustrates the joining of two unconnected endpoints using a shortest-path algorithm. Figure 2D illustrates resolution of an unconnected endpoint, using extrapolation of the endpoint's trajectory. Full details of the skeleton correction methods can be found in the Supplementary Methods. 

Figure 2

 

A suite of correction methods, both automated and manual, were employed to achieve accurate segmentation (see Supplementary Methods for full details). (A) The original skeleton (red) is overlaid on the binary mask (gray), with corrections (determined using the methods in [B–D]) in green. (B) Skeleton correction by automatically combining skeletons at different thresholds. The bridging strand (green) is not present in the initial skeleton (red), but at a lower threshold level the binary mask is extended by a region of lower intensity, which after skeletonization results in the connection. (C) Skeleton correction by joining endpoints via a shortest-path computation. The manually identified endpoints are joined together automatically using a method based on Dijkstra's algorithm,³⁹ where weights, computed from the enhanced image data, are used to determine the optimal bridging strand. (D) Skeleton correction by extrapolation of an endpoint's trajectory. The bridging strand is extrapolated from the start point, based on the local curvature and intensity, until it reconnects with the skeleton.

Vascular Graphing

To construct a graph of the vasculature, the corrected skeleton was first converted to a vectorized map. Vectorization entails explicitly labeling all branches and nodes (branch points) that constitute the corrected skeleton, as well as determining the interconnectivity between them (see Supplementary Methods). We described the branching hierarchy of the retinal blood vessels using branch order, which describes the minimum number of bifurcations to reach a given vessel branch, assuming an order of “1” for major vitreal vessels. Vessel branch length and diameter were estimated using a 3D parametric curve-fitting approach (see Supplementary Methods for full details). 

Flow Simulation for Arteriole and Capillary Dilation

A model-based reconstruction of blood flow, using a resistive network, was used to further probe the interconnectivity of the retinal vasculature (see Supplementary Methods). Absolute values of flow and pressure are not necessarily accurate in such a model, as they depend on assumptions relating to the boundary conditions and blood viscosity. However, the distributions of these parameters are largely independent of such assumptions,⁴⁰ and thus flow simulations based on the model can be used as a means to conceptualize vascular interconnectivity. 

In order to investigate how the topology and connectivity of the vascular network impact blood flow control, dilations of selected arterioles, and capillaries were simulated. The purpose of these simulations was to investigate whether local retinal blood flow could plausibly be controlled through localized vasodilation. Dilation was simulated by increasing the diameter of a particular vessel by an amount on the order of experimentally measured dilations in rats.^2,3 Two types of dilation experiments were performed. Firstly, simultaneous dilation of all vessels of a particular category (either arterioles or capillaries) was performed in order to elicit a bulk response in flow. For arterioles, the specified diameter increase was 10% to 12%, whereas for capillaries it was 0% to 2.5%.^2,41 The proportionate change in total capillary flow across each vessel layer, in response to dilation, was computed. Secondly, a single arteriole was dilated by 10%. This procedure was repeated for every arteriole in the network, and the average proportionate change in flow as a function of branching distance from the dilated arteriole was computed. 

These simulations presented a unique opportunity to assess retinal blood flow regulation; however, several limitations should be mentioned. The localized dilation of a single arteriole may not be realistic, and beyond the limits of the spatial precision of neurovascular coupling in the retina.² Also, our assumption of constant flow boundary conditions (see Supplementary Methods) may not be realistic. While the arteriole to dilate was chosen away from the boundary of the imaged region, the assumption of constant flow boundary conditions may have required more redistribution of flow than would have been necessary during realistic functional hyperemia in the retina. In reality, the flow supplying and draining the retinal sector with arteriolar dilation would probably increase, thus leading to a more uniform increase in flow than was predicted by our simulations (which, because the boundary conditions were not changed, implicitly required balancing the local increase in flow with a decrease in flow elsewhere). 

Results

Topologic features of the retinal vasculature were derived from the skeletonized data (Fig. 3). We corrected and skeletonized the inner retinal vasculature, and delineated each of the vascular layers by color (Fig. 3A; see also Supplementary Movie S1). The proportion of total skeleton branches that were introduced by the skeleton correction procedures (Fig. 2) was 18.4 ± 4.0% (N = 3). We then determined the averaged linear vessel density in each layer, computed as the total length of all vessels in the layer divided by the field of view size (Fig. 3B), and compared the distribution of vessel diameters by layer (Fig. 3C). A Kruskal-Wallis test⁴² was used to determine if the branch diameters in each of the three different layers were similarly distributed. The test rejected the null hypothesis that the branch diameters in each layer had the same distribution, at a 1% significance level. Thus, vascular metrics differ in each layer, pointing to the varied biology of inner retinal vasculature. 

Figure 3

 

Topology of the retinal vasculature. (A) Skeletonization of inner retinal vasculature with different layers denoted by color (see also Supplementary Movie S1). (B) Linear vessel density for the three vascular layers. The error bars represent the SEM. (C) Histogram of estimated vessel branch diameters, computed over three separate retinal datasets. A Kruskal-Wallis test rejected the null hypothesis that branch diameters were similarly distributed for the three vascular layers (P < 0.01).

The graph representation further enabled laminar analysis of branching order (Fig. 4). We determined the order of all vessel branches in the graph (Fig. 4A, see also Supplementary Movie S2), by assigning an order of 1 to the major supplying and draining vessels and then using the connectivity information from the graph to infer the order of all other branches (as illustrated in Fig. 4B). We then compared the distribution of branch order among the three vascular layers (Fig. 4C). A Kruskal-Wallis test was used to determine if the branch orders in the three different layers were similarly distributed. The test rejected the null hypothesis that the branch orders in each layer had the same distribution, at a 1% significance level. Thus, inner retinal microvasculature, supplied and drained through the vitreal arteries and veins, respectively, forms a “tree,” with higher order branches in deeper layers. 

Figure 4

 

Laminar analysis of branching orders. (A) Map of vessel branch order (see also Supplementary Movie S2); (B) diagram illustrating the definition of branch order: major supplying arteries and veins have an order of 1, arterioles and venules have an order of 2, and other vessels (capillaries) have an order of 3 or higher; (C) Histogram of vessel branch orders, illustrating the distribution of branches of different order among the three vascular layers. A Kruskal-Wallis test rejected the null hypothesis that branch orders were similarly distributed for the three vascular layers (P < 0.01).

Laminar analysis of interconnectivity of the 3D vascular network was performed through highlighting connecting vessels (Fig. 5). The density of vessels connecting between the superficial and intermediate layers for the illustrated dataset was 77.2 mm⁻² while the density of vessels connecting between the superficial and deep layers was 26.9 mm⁻². In contrast, the vasculature in the intermediate and deep layers (the plexiform layers) was intimately connected, with an interconnecting branch density of 253.0 mm⁻². Supplementary Figure S1 shows the density of offshoots (arterioles and venules) per unit length, for the major supplying/draining arteries and veins. Implications of the asymmetric connectivity for blood flow regulation are investigated below. 

Figure 5

 

Laminar vascular interconnectivity. The maps highlight (red) vessel branches whose two associated node points are located in different vascular layers (i.e., branches that connect different vascular layers to each other).

While the graph is informative, its utility for inferring blood flow control mechanisms is limited. Therefore, we further developed a simple network flow model, informed by experimental boundary conditions, based on the graph. Initialization of the flow simulation was based on the measured Doppler OCT flow rates, assigned to the large vessel branches in the superficial layer (see Supplementary Methods). Experimental flow rates were not available for loose vessel branches in the capillary layer. In the absence of other information, these branches were assigned a flow rate ascertained by averaging the absolute flow rates computed for capillaries of similar diameter and assigning the flow direction randomly. 

We then mapped the simulated effects of vessel dilation, by computing proportionate changes in flow for each vessel (Fig. 6). Figure 6A shows a map of the normalized estimated change in flow resulting from simulated dilation of an arteriole branch. The dilation was simulated by increasing the arteriole's diameter by 10% and recomputing the flow in the network. Figure 6B shows an oblique 3D view of the flow change map. This result confirms the capacity of the inner retinal vascular network for local blood flow regulation. 

Figure 6

 

Spatial localization of changes in flow after simulated vessel dilation. (A) Map of the normalized change in flow produced by a 10% dilation of an arteriole branch. The dilated branch is marked in green; (B) zoomed-in oblique 3D view of the map in (A).

Figure 7 summarizes both the laminar variation and spatial localization of flow changes after simulated vessel dilation. We investigated dilation of capillaries and arterioles, both candidates for active control of blood flow. We first computed proportionate changes in total capillary flow in the three retinal layers when all vessels in a given category were simultaneously dilated by random amounts based on previously reported experimental measurements.² Five repeated simulation experiments were performed. The averaged results and standard errors are shown in Figure 7A. We then computed the proportionate change in total flow due to simulated dilation of an individual arteriole by 10%. The dilation experiments were repeated for all arterioles in the network, choosing a different arteriole each time. The averaged results and standard errors are shown in Figure 7B, as a function of branch order relative to the dilation site (as illustrated in Fig. 7C). In summary, arteriole dilation causes localized and lamina-specific flow changes, while bulk dilation of capillaries in a given layer effectively recruits flow to that particular layer. 

Figure 7

 

Laminar dependence and spatial localization of flow changes after simulated vessel dilation for various vessel categories: arterioles (superficial-deep; ASD), arterioles (superficial-intermediate; ASI), capillaries (superficial; CS), capillaries (intermediate; CI), capillaries (deep; CD), capillaries (connecting intermediate to deep; CC). (A) Proportionate change in the total flow among all capillaries in the three retinal layers, due to bulk dilation of a given vessel category; (B) response to a single 10% arteriolar dilation, as a function of branch order relative to the dilation site and averaged over all arteriolar dilation sites; (C) designation of branch order relative to dilation site. The dilated vessel in the illustrated case is an arteriole connecting the superficial and intermediate layers.

Discussion

Our results demonstrate, to our knowledge, the first 3D vectorized graph of the inner retinal vasculature. Our results confirm that the inner retina is supplied from the large vitreal arteries to a “trilaminar” hierarchy of higher-order vessels consisting of the ganglion cell layer/nerve fiber layer (superficial), inner plexiform layer (intermediate), and outer plexiform layer (deep). The intermediate and deep layers of vasculature were found to stratify with dendrites in the plexiform layers, where metabolism is highest.⁴³ Particularly, the intermediate layer was located around the off-sublamina at the outer edge of the inner plexiform layer.⁴⁴ The three layers are drained via the vitreal veins. 

The graphing methodology demonstrated in this work has provided insight into the inner retinal architecture beyond the well-accepted description above. In particular, a branching tree parallels the laminar hierarchy, with lower order vessels on average in the superficial layer, compared with higher order vessels on average in the deep layer, with the intermediate layer tending to have more intermediate branching order vessels (Fig. 4C). Moreover, there are clear laminar differences in vessel density, with the deep layers having the highest density (Fig. 3B). In view of the fact that the deep layer also has the highest vessel orders (Fig. 4C), the higher density in the deep layer may be viewed as compensating for lower flow per vessel and the presumed lower oxygenation in higher order vessels. 

The capacity for laminar control of blood flow in the retina depends critically on the nature of connections between the inner retinal layers, which were analyzed for the first time in this work. The interconnectivity between the superficial and plexiform layers was found to be relatively sparse, whereas the intermediate and deep capillary layers have an abundance of interconnecting vessels (Fig. 5). The sparse connectivity between the superficial vasculature and the capillary plexuses may carry significance for disease, as the plexiform layers are among the most oxygen-consuming layers in the retina.¹ The implications of this pattern of inner retinal vascular connectivity on susceptibility to ischemia will be the subject of future investigations. 

Flow simulations based on the graph further demonstrated the capacity of the inner retinal vascular architecture for both laminar and localized blood flow control. Simulated results based on a network flow model showed that flow changes in a particular capillary plexus can be achieved by dilation of its feeding arterioles (Fig. 7A). Bulk dilation of capillaries in a given plexus also preferentially recruited flow to that plexus. The redistribution of flow resulting from a dilation is local, with changes diminishing sharply away from the dilation site (Figs. 6B, 7B). In some cases, decreases as well as increases in flow were observed in the vessels surrounding the dilation site, ostensibly due to a passive redistribution of flow around the activated vessel (Figs. 6B, 7B). The spatial flow response patterns observed during simulated dilation are largely consistent with the interconnectivity of the network. 

Though flow simulations were useful to assess connectivity, a number of technical limitations should be mentioned. Errors in baseline branch diameters can cause large errors in simulated flow. To improve diameter estimates, higher transverse resolutions will be enabled with adaptive optics in the future. Arterioles (defined here as first-order branches from vitreal arteries) and capillaries were singled out for simulated dilation, as they are known or hypothesized to be involved in blood flow regulation.⁴⁵ The dilation was modeled as a fixed increase, not localized to any particular part of the vessel branch, and temporal effects were not considered. Results are also qualified by insufficient knowledge of how boundary conditions may change under dilation. With these caveats in mind, the flow simulations confirmed the capability of the trilaminar architecture to perform both laminar and spatially localized control of inner retinal blood flow. More complete physical models of flow will be investigated in the future. 

Conclusions

The model of the rodent inner retinal vasculature presented in this work is, to our knowledge, unique in its level of topologic detail. Vascular branching patterns within the capillary plexuses are represented along with supply and drainage routes. The capability for laminar and localized control of blood flow is theoretically confirmed by treating the vasculature as a resistive network. Models that integrate vascular structure and flow distributions in 3D will advance understanding of the relationship between blood flow, metabolism, and neuronal activity, and provide a framework for characterizing the progression of retinal diseases. 

Acknowledgments

Supported by the National Institutes of Health (R00-NS067050, R01NS094681, P30AG010129; Bethesda, MD, USA), American Heart Association (11IRG5440002; Dallas, TX, USA), the Glaucoma Research Foundation Catalyst for a Cure Biomarker Initiative (San Francisco, CA, USA), and Research to Prevent Blindness, Inc. (New York, NY, USA). 

Disclosure: C. Leahy, None; H. Radhakrishnan, None; G. Weiner, None; J.L. Goldberg, None; V.J. Srinivasan, P