A Novel and Open-Source Application for Vasculature Dataset Analysis and Visualization

Vascular networks can dictate and indicate states of health and disease. Structural analyses of 35 these networks can facilitate improved understanding of disease states. Recent advances in preclinical 36 imaging techniques and segmentation software have led to the generation of large-scale vasculature 37 datasets. However, these advances have not been accompanied by the development of modernized, 38 open-source analysis software packages. Here we describe VesselVio, an application developed to 39 analyze and visualize pre-segmented 2D and 3D vasculature datasets. Vasculature datasets can be 40 loaded and analyzed with custom parameters to extract numerous quantitative whole-network and 41 individual segment features. Visualization of results and accuracy inspections can be conducted using 42 the interactive visualization tool. The utility and compatibility of VesselVio is demonstrated via the 43 analysis of 3D inferior colliculus segmentations from female and male mice as well as the analysis of 44 2D retinal fundus images of control, glaucomatous, and diabetic retinopathy patients. 45 46 47 48 49


Centerline Extraction and Graph Theory Enable Detailed Vascular Network Feature Extraction 91
We sought to build an open-source application that allows users to extract numerous quantitative 92 features from 2D and 3D vascular networks. Different imaging techniques can produce varying qualities 93 of raw images that require unique vessel detection and thresholding processes. Because of this, 94 VesselVio was developed for vascular datasets that have already been binarized and segmented, 95 enabling it to be used with datasets of any imaging origin ( Figure 1A). This development choice enabled 96 analysis pipeline optimization for increased accuracy and speed of feature identification. The 97 visualization component for the program enables researchers to examine their results alongside original 98 volume meshes for accuracy inspection. 99 The features of a vascular network can be extracted by identifying centerlines of the network 00 and creating undirected graphs of the centerline points 23,24 . To locate vessel centerlines, a widely used 01 2D and 3D thinning algorithm for binarized images was employed to extract the skeleton of our 02 datasets 25 . Edge-connectivity of the skeleton points is identified and used to create an undirected graph 03 with the resulting edges and points, G = (V, E). Following initial graph construction, the datasets follow 04 a series of correction processes to filter out spurious branchpoint labels, smooth centerlines, and 05 remove isolated/endpoint segments at user-defined size ( Figure 1B). These initial processing stages 06 enable downstream quantifications of network and segment features ( Figure 1C). 07 datasets are loaded into the program where analyses can be conducted using custom analysis and 09 feature export parameters. b, Centerlines are extracted from the loaded datasets, and undirected 10 graphs are created with the resulting vertices and edges. Centerlines are smoothed, spurious branch 11 points are filtered, and segments are filtered prior to result extraction. c, The resulting datasets can 12 be visualized for inspection of the accuracy and quality of results. 13 14

Endpoints and Branchpoint Clique Filtering 16
After graph construction, network endpoints and branchpoints are identified by examining the 17 degree of connectivity of the centerline vertices ( Figure 2A). Centerline neighbor 26-connectivity is a 18 commonly used method for branchpoint and endpoint identification 7,16,18 . Although endpoints can be 19 accurately detected using this technique, branchpoints are spuriously labeled and counts become 20 artificially inflated ( Figure 2B). Some programs allow for interactive user input to correct incorrectly 21 labelled segments, centerlines, or branchpoints (e.g., Imaris 8 ). However, manually correcting 22 branchpoints in 3D datasets can become tedious and time-intensive, particularly when datasets are 23 several gigabytes large. Because VesselVio is designed for automated analysis without user guidance, 24 a set of algorithms were implemented to automatically filter incorrectly identified branchpoints from 25 analyzed datasets. 26 Spuriously labelled branchpoints form small loops in the constructed graphs, otherwise known 27 as cliques. Branchpoint cliques can be identified by filtering the graph of endpoints, 2-degree vertices, 28 and single or two-component branchpoints. These cliques contain candidate branchpoints with one or 29 more loops generated by between-edges ( Figure 2A) method. *** -p < 0.001, **** -p < 0.0001. 56

Identifying Segment Radii Based on Voxel Edges 58
Several techniques exist for identifying vessel radii. One technique involves recording the largest 59 maximally inscribed spheres that can rest within mesh vessel centerline points 29,30 , but this technique 60 often depends on the creation of directed graphs (i.e., manually directed vessel hierarchy) and thus 61 was not suitable for an automated pipeline. A similar method identifies the Euclidian distance between 62 a vessel centerline and the center of the nearest non-vessel neighbor 6,31 . However, one pitfall to this 63 approach is that vessels with near-resolution or at-resolution radii are incorrectly sized when their 64 closest non-vessel neighbor is located along the 1 st dimension of individual axes ( Figure 2C), leading 65 to oversized single-pixel/-voxel vessels ( Figure 2D). This issue is not as apparent for 2 nd and 3 rd 66 dimension neighbors (Supplementary Figure 1). As such, a simple half-unit correction for 1D neighbors 67 (6-connectivity) is implemented to preserve small-segment radii measurements ( Figure 2C). This To demonstrate VesselVio 3D compatibility, several datasets were analyzed. First, 85 cerebrovascular resin casts of female and male CFW mice (n = 5) were created as described 86 previously 8 . Following tissue removal and subsequent osmication, the resulting vasculature casts were 87 imaged using a µCT scanner at an isotropic resolution of 2.7 µm 3 . Next, the inferior colliculi (IC) were 88 segmented and binarized as example datasets (Supplementary Video 1). As expected, the IC were 89 densely vascularized with large peaks in small diameter vessels ( Figure  length (e), mean segment tortuosity (f), and mean segment length (g). h, Distribution of the number of 08 segments along 0-20+ radius bins. i-j, example images from female (i) and male (j) inferior colliculus 09 reconstructions. Results for a-d and h were corrected by the region-of-interest volume for individual 10 segmentations; these results are per mm 3 . Data are represented as mean ± SEM. a-g, Data analyzed 11 using two-tailed student's t-test. h, Data analyzed using 2-way ANOVA; n = 5. ** -p < 0.01. 12 13 2D Retinography 14 To demonstrate the utility of VesselVio with 2D datasets, retinographs sourced from the HRF 15 Image Database were analyzed and compared 33 . This database contains high-resolution images of 16 healthy control patients, patients with diabetic retinopathy, and patients with glaucomatous eyes.    The first set of analyses examined sex-differences in the cerebrovasculature of female and male 51 CFW mice. There are numerous known sex-differences in the cardiovascular system, including 52 vascular tone 35 , microcirculation 36 , and blood-brain barrier function 37,38 . We sought to examine how the 53 cerebrovascular network in a specific nucleus, the inferior colliculus (IC), may differ based on sex. 54 Analysis revealed few overall differences in the vasculature of these nuclei, save for differences in 55 average segment lengths. However, these results still serve to demonstrate 3D network analysis utility 56 of VesselVio. 57 The second set of analyses examined retinal vascular networks of healthy controls, DR patients, 58 and patients with glaucoma. Alterations in retinal microvascular networks are associated with vision 59 loss 34 , and structural changes associated with excessive or reduced angiogenesis can also serve as 60 indicators for underlying disease states 39,40 . Our analyses revealed numerous vascular differences 61 among the groups, including alterations in branchpoint density, altered vascular area, and altered 62 tortuosity in the two disease states. Our analyses recapitulated previously reported reduced vessel 63 area and diameter in glaucomatous eyes 41 , as well as increased segment tortuosity of DR eyes 42 . 64 Contrasting results in DR eyes were noted compared to previous studies that reported increased vessel 65 diameters 43,44 , whereas we and others observed decreased diameters 45,46 . However, these differences 66 may be due to the vessel type (the HRF vessels are not separated by arterial/venous hierarchy) or 67 vascular region analyzed. Although VesselVio is not intended to be used in any clinical diagnostic 68 contexts, this analysis demonstrates the ability of this application to identify and characterize vascular 69 network alterations in pathological conditions. 70 In common with most analytical program, the functionality of VesselVio is limited by the quality 71 or resolution of the images that are loaded into the program. For example, if anisotropic datasets are 72 loaded without pre-smoothing or blurring filters, then skeletonization of these datasets can produce 73 erroneous segments. To adjust for these errors, the option to prune small, connected end-segments is after perfusion, the craniums were decalcified with a 12-h wash of 5% formic acid (BDH4554; VWR 08 International, USA), the brains were dissected, and remaining tissue was removed from the casts with 09 two 12-h washes of 7.5% KOH (BDH7622; VWR International) at 35 ºC. Casts were then rinsed with 10 three 1-hour Milli-Q water washes, and the cleaned casts were osmicated in a 1% solution of osmium 11 tetroxide (#75632; Sigma Aldrich, USA) for 12-hours to allow for optimal x-ray diffraction during µCT 12 scans. Casts were imaged on a SkyScan 1272 (Bruker, USA) at 50 kV/ 200 µA with 360º rotations in 13 step sizes of 0.17º, no filter, 900 ms frame exposures, and 4 frame averages/step to produce an 14 isotropic voxel resolution of 2.7 µm 3 . Scan parameters were determined based on the manufacturers 15 guidance to achieve optimal x-ray transmission through the sample. Scans were then reconstructed 16 using NRecon (Bruker) with beam hardening corrections at 15%, ring artefacts reduction at 3,

Dataset Input and Processing Preparation 25
Images are loaded into VesselVio using the Simple-ITK image reader to improve compatibility 26 with various file formats 48 . It is important to note that VesselVio is only compatible with vascular datasets 27 that have been pre-segmented and binarized. To ensure that datasets loaded into the program are 28 prepared appropriately for subsequent analysis, we re-binarize all inputs with a threshold of 1 and 29 convert the image stack or single image into a binarized contiguous array of 0-value background and 30 1-value foreground elements. If images are not loaded in 8-bit grayscale, then this binarization to 31 unsigned integers process can save memory during subsequent processing. All array processing in 32 VesselVio is conducted with NumPy 49 . 33

Volume Skeletonization and Centerline Extraction 35
VesselVio employs the scikit-image 50 implementation of a widely used medial axis parallel 36 thinning algorithm 25 . We selected this algorithm because it produces few erroneous branchpoint 37 extensions, particularly when used with high-resolution datasets. This algorithm is also capable of 38 thinning 2D and 3D datasets, making it optimal for analytical purposes. Following skeletonization, (n,  Figure 1). However, the traditional EDT method will record this vessel as having a 52 radius of 1 µm 3 , doubling the apparent diameter. This is because ED measurements are calculated 53 between the coordinates of the center of the centerline voxel and the center of its non-vessel neighbor. 54 Radii inflations in traditional EDT measurements are less pronounced as vessels become larger, but 55 small-diameter vessels can be binned incorrectly using this blanketed technique. Instead, the 56 calculation should be made to the edge of the nearest non-voxel neighbor if it is oriented single-57 dimensionally to the centerline point. Because of this, a modified EDT (mEDT) calculation was 58 implemented that has 0.5-unit corrections for all neighbors oriented along single X, Y, or Z axis 59

directions. 60
An argument could be made to calculate distance to the corner of the center pixel/voxel, but the 61 edge was selected. Given the phenotypical tubular shape of vessels, EDT-based radii measurements 62 are more appropriate for neighbors located along two-and three-dimensional orientations in 2D or 3D 63 space. Thus, only radii measurements for nearest neighbors in 1D space are corrected (Supplementary 64 identification process prevents parallel edges from being created 6 . 98 (1) (2) Following graph creation, n!=2 degree vertices are filtered from the graph. The remaining 99 components are then scanned, and isolated segments and endpoint segments shorter than user-00 defined sizes are pruned (Supplementary Video 5). Following the removal of isolated segments from 01 the graph, a custom flood-filling process is implemented to remove corresponding pixels/voxels from 02 the original volume. Flood-filling of these isolated segments prevents inflation of the calculated 03 area/volume. This technique is not implemented for endpoint segments, as endpoint segment pruning 04 is used to remove erroneous centerlines created during skeletonization. From our tests, endpoint 05 segments in need of pruning are often incorrectly skeletonized centerlines within large segments, so 06 flood filling is not implemented for these segments. 07 Graphical representation of the vasculature construction also allows for improved identification 08 of branchpoints. Branchpoints are often identified based on simple neighborhood connectivity; if a 09 centerline point has more than two neighbors, it is defined as a branchpoint 7,16 . However, this method 10 leads to an artificial inflation of identified branchpoints; this is because at branchpoint junctions, multiple 11 vertices can have > 2 neighbors -not just a single branchpoint. In graphical space, these falsely 12 identified branchpoints form small loops, or cliques. To correct for this, we constructed a set of 12 13 branchpoint filtering algorithms that remove >99 % of identified branchpoint cliques as determined by

Mesh Visualization and Application Interface 44
To construct meshes for visualization, we leverage the high-level python package PyVista 56 that 45 wraps The Visualization Toolkit. We create individual poly-datasets from our segment splines, apply 46 tube filters to create simple network and scaled network (based on average radius) meshes, and assign 47 each tube a radius, length, and tortuosity scalar for visualization. All scalars and scaled segment sizes 48 are based on the mean of segment features. These segments are combined into an undirected grid for 49 surface extraction and subsequent rendering. Additional branchpoint and endpoint meshes are created.  Figure 2). 55 56

Statistical Analysis 57
All statistical analyses were conducted using Prism 9 (GraphPad; USA). Synthetic vasculature 58 results were analyzed using one-way ANOVA tests and repeated-measures two-way ANOVA tests. 59 HRF datasets were analyzed using one-way ANOVA tests. Mice IC data were analyzed using two-60 tailed student's t-tests. Distributions of segment counts per radii bin, average segment length per radii 61 bin, and segment tortuosity per radii bin were all analyzed using two-way ANOVA tests. Following main 62 effects observations in one-or two-way ANOVA tests, multiple comparisons were made using Tukey's 63 test. An alpha of 0.05 was set for statistical significance.