Micrometer-resolution reconstruction and analysis of whole mouse brain vasculature by synchrotron-based phase-contrast tomographic microscopy

Nervous tissue metabolism is mainly supported by the dense thread of blood vessels which mainly provides fast supplies of oxygen and glucose. Recently, the supplying role of the brain vascular system has been examined in major neurological conditions such as the Alzheimer’s and Parkinson’s diseases. However, to date, fast and reliable methods for the fine level microstructural extraction of whole brain vascular systems are still unavailable. We present a methodological framework suitable for reconstruction of the whole mouse brain cerebral microvasculature by X-ray tomography with the unprecedented pixel size of 0.65 μm. Our measurements suggest that the resolving power of the technique is better than in many previous studies, and therefore it allows for a refinement of current measurements of blood vessel properties. Relevant insights emerged from analyses characterizing the regional morphology and topology of blood vessels. Specifically, vascular diameter and density appeared non-homogeneously distributed among the brain regions suggesting preferential sites for high-demanding metabolic requirements. Also, topological features such as the vessel branching points were non-uniformly distributed among the brain districts indicating that specific architectural schemes are required to serve the distinct functional specialization of the nervous tissue. In conclusion, here we propose a combination of experimental and computational method for efficient and fast investigations of the vascular system of entire organs with submicrometric precision.


Abstract 24
Nervous tissue metabolism is mainly supported by the dense thread of blood vessels which mainly provides 25 fast supplies of oxygen and glucose. Recently, the supplying role of the brain vascular system has been 26 examined in major neurological conditions such as the Alzheimer's and Parkinson's diseases. However, to date, 27 fast and reliable methods for the fine level microstructural extraction of whole brain vascular systems are still 28 unavailable. We present a methodological framework suitable for reconstruction of the whole mouse brain 29 cerebral microvasculature by X-ray tomography with the unprecedented pixel size of 0.65 μm. Our 30 measurements suggest that the resolving power of the technique is better than in many previous studies, and 31 therefore it allows for a refinement of current measurements of blood vessel properties. Relevant insights 32 emerged from analyses characterizing the regional morphology and topology of blood vessels. Specifically, 33 vascular diameter and density appeared non-homogeneously distributed among the brain regions suggesting 34 preferential sites for high-demanding metabolic requirements. Also, topological features such as the vessel 35 branching points were non-uniformly distributed among the brain districts indicating that specific architectural 36 schemes are required to serve the distinct functional specialization of the nervous tissue. In conclusion, here 37 we propose a combination of experimental and computational method for efficient and fast investigations of 38 the vascular system of entire organs with submicrometric precision. 39 40

Introduction 41
Cerebral blood vessels sustain neuronal activity by providing metabolic components and oxygen to the brain 42 tissues and by removing catabolic waste products. Specifically, it has been estimated that in humans and 43 primates, synaptic activity and action potentials account for about 96% of the total energy consumed 1 , events 44 enabled by a tight coupling among neuronal, glial and vascular cells. While much efforts have been spent to 45 reconstruct in detail the sophisticated networks of specific neuronal circuits, no comparable achievements have 46 solved the intricate microvascular architecture even in small nervous systems. Precise depiction of the vascular 47 structure is also important for the comprehension of the biophysical mechanisms governing the crucial 48 interplay among neurons, glial cells and blood vessels, a physiological event known as neuro-vascular 49 coupling 2 . 50 The brain vascular structure is altered in many neurological diseases such as cerebrovascular diseases and 51 various forms of Alzheimer's and Parkinson's diseases which constitute major health issues worldwide 3 . Many 52 of these conditions are related to changes in the structure of the cerebral blood vessels, such as increased or 53 decreased capillary density, microbleeds, stiffening of artery walls, or increase in vessel tortuosity 4,5 . 54 Reconstructing the structure of the whole vascular network in a brain is challenging as the whole brain is 55 generally very large compared to the smallest microvessels. Two imaging techniques often applied to this 56 purpose, in the context of mouse brain analyses, are magnetic resonance imaging 6 and light sheet microscopy 7-57 9 . In these techniques resolution and pixel size are however often closer to 10 μm rather than 1 μm, and usually 58 too large to unveil the smallest microvessels that have diameters down to less than 5 μm. On the other hand, 59 X-ray tomographic microscopy (CT) offers higher resolution and the advantage of not requiring extensive 60 sample preparation. Preparations such as tissue clearing increase the possibility of deformations in the sample, 61 and their efficacy could vary spatially 10,11 . 62 Previously the cost of high resolution in CT was small sample size 12 . However, recent advances in imaging 63 speed 13,14 and post-processing algorithms 15 of these are not openly available to the  70  whole scientific community, are in early stage of development, or are geared towards single technique or  71 imaging modality. 72 In this work, we show that synchrotron-based phase-contrast tomographic microscopy can be advantageously 73 used to image micro-vessels in the whole mouse brain. In our technique the microvessels were perfused with 74 contrast agent (Indian Ink 19 ) and fixed with formalin solution. The perfused sample was tomographically 75 investigated in a mosaic imaging mode, and the individual tiles were combined into one large volume image 76 of the whole brain using a non-rigid stitching algorithm 15 . A main advantage of the mosaic imaging mode is 77 that slow deformations of the sample do not lead to imaging artifacts, in contrast with more traditional CT 78 imaging techniques. Finally, the full 11 TB volume image was segmented and analyzed using an in-house 79 developed and freely available software, capable of processing the images in a few days using a small computer 80 cluster as detailed below in the Section 4. 81 From the volume images we generated a vascular graph consisting of vessel branches and bifurcation points. 82 The graph embedded various information related to the topology and morphology of the microvessel branches, 83 e.g. length and average diameter. Subsequently, we determined the anatomical regions corresponding to our 84 blood vessel space by registration of the volume image to the Allen mouse brain atlas 20 . Eventually, we 85 described the vascular anatomy of the brain in different anatomical regions using quantities such as 86 microvessel tortuosity, bifurcation density, vessel length density, and intercapillary distance. 87 88 2. Results 89

Sample preparation and imaging 90
Whole-brain samples were prepared by intracardiac perfusion of the vascular system with Indian Ink, 91 according to the protocol in 19,21 . The brain was then extracted and stored in polyphosphate buffered solution 92 to maintain constant hydration until and during the data acquisition. The brain sample was CT imaged in a 93 mosaic phase-contrast imaging mode where partially overlapping individual tomograms are taken side-by side 94 to cover the whole volume of the brain. The tomograms were stitched into one large volume image using a 95 non-rigid stitching algorithm 15 that accounts for small deformations of the overlapping regions. The 96 dimensions of the final volume image were approximately 13600 × 13600 × 30000 pixels with a pixel size of 97 0.65 µm. Finally, the volume image was denoised and segmented using standard image analysis algorithms. 98

Image quality 99
Initially, the quality of the volume image and the segmentation was visually evaluated by several imaging 100 specialists. The algorithmic segmentation was found to visually match most of vessels identifiable on the 101 volume image, see Figure 1 and  Figure 1A-D 107 shows various visualizations of the segmented microvessels. 108 Based on visual inspection of the blocks that were identified as badly segmented, the largest segmentation 109 errors were found near large arteries and veins (diameters typically in the range of tens or hundreds of 110 micrometers) that are not thoroughly perfused with the contrast agent (Indian Ink). We have not considered 111 these vessels in the analysis. However, as shown by the statistics calculated from the confusion table, the 112 number of such vessels is small compared to the microvessels, and therefore they do not contribute 113 significantly to the results. 114

Overall structure of microvessels 115
In order to quantify the structure of individual microvessel branches, we found the centerline of each branch, 116 and the bifurcation points of two or more branches. Additionally, we measured the length and the diameter of 117 each branch. 118 We estimated the whole brain to contain approximately 4. bifurcation might lead to varying estimates of bifurcation density. Length density did not suffer from such an 127 ambiguity in its definition and its value agrees better with the previously reported ones. 128 For the first time, our results returned estimations of the whole-brain vascular length (295 m), and of the 129 average vessel branch length (53 ± 3 μm) ( Figure 1G) at the level of the single microvessels. Average vessel 130 diameter was 5.8 ± 0.4 μm in the whole brain ( Figure 1H). It is conveniently between previously reported 131 values of 4.25 μm 22 and 8 μm 7 . The differences might be caused by various sample preparation routines such 132 as perfusion and optical clearing. It was not certain that possible shrinkage or swelling of the vessels caused 133 by these operations can be easily accounted for 24 , particularly when the focus is on local micro-scale properties 134 and not on overall average deformation. In particular, any inaccuracy caused by local non-isotropic 135 deformations are easily propagated to the results due to the small diameter of the vessels (in pixels, 5.8 μm = 136 8.9 pixels). According to results in Figure 1L vessel diameter was almost constant in the longer branches but 137 varied more in the shorter ones. 138 The length over diameter ratio (L/d, dimensionless length) is related to the pressure loss in the vessel through 139 the Darcy-Weisbach equation. Large values of L/d indicate that the vessel is long and thin, and the flow 140 resistance of the vessel is large. We found an average L/d value of 11 ± 1 ( Figure 1I) and largest values ranging 141 to more than 100. The values indicate a large spectrum of putative flow resistance regimes 25 . 142 Tortuosity is a measure of how much a blood vessel segment twists, with high values typically related to 143 pathologies 26 . Average tortuosity of the vessels was 1.24 ± 0.01 ( Figure 1J), a value which was maximized 144 with 20-30 µm length vessels and generally increases with the vessel length ( Figure 1M). Further, tortuosity 145 was highest in vessels of diameters between 5 to 7 µm ( Figure 1N), indicating that the purpose of the vessels 146 in this size range is to assist in even transport of metabolic components and waste products to and from the 147 tissue, in contrast to efficiently transferring them for long distances. The average distance the products must 148 transport outside of blood vessels equals to the distance to the nearest microvessel, and that was measured to 149 be 15 ± 1 μm. The value corresponds to approximately 2.5 average microvessel diameters ( Figure 1K). 150

Differences in the vascular structure between anatomical regions 151
In order to obtain a more detailed picture of the structure of the vessel network, we co-registered the obtained 152 mouse brain vessels with the Allen mouse brain atlas at the highest resolution of the atlas (100 μm pixel size). 153 In addition, we clustered the atlas brain regions into two different hierarchical groups such that very small 154 regions were combined in order to guarantee that each clustered region contained statistically significant 155 number of vessel branches. In the first (finer) level, we have 44 different regions (Supplementary table 1, left  156 column, regions), and in the second coarser level 11 regions (Supplementary table 1 Figure S17). In conclusion, geometrical, 187 morphological and topological features of mouse brain microvessels appeared regionally specific suggesting 188 distinct roles in support of local specialization of brain districts. 189

Discussion 190
The results highlighted peculiar characteristics of specific macroregions that were mainly the white matter, the 191 olfactory bulb and the cerebellum but also the striatum and the somatosensory, motor and visual cortices, 192 which appeared to get extreme values in our estimations. Indeed, we observed a strong correlation between 193 the neuronal density 28,29 and the numbers of branch points and tortuosity ( Figures 4A and K), a weaker but 194 sustained correlation has been detected also with the distance to the nearest vessel ( Figure 4N). In addition, 195 comparable and in some cases even stronger correlations held also in the glial cell density (Figure 4, third and  196 fourth columns). Note that the estimations were biased by the zero neuronal density of white matter and, 197 oppositely, by the high neuronal density of cerebellar layers (more than fivefold the average of other regions). 198 Indeed, removing these two regions from the linear regressions results in stronger correlations between 199 neuronal density and many of the measured quantities (Supplementary figure 2). Such postliminary 200 considerations suggest that no trivial rules govern microvascular features in relation to the other existing 201 cellular families (neurons and glia). Surprisingly, vessel tortuosity appeared to be the best predictor (as for 202 linear regression) for neuronal density while the distance to the nearest vessels played the same role for glial 203 density. 204 Besides density correlations, a reasonable observation was that typical high energy demanding regions were 205 denser of vessels with long segments and small diameters. A set of neurophysiological and anatomical 206 considerations supported the observed results. The brain white matter is mostly composed by myelinated 207 neuronal axons and glial cells. From metabolic perspective, axonal segments demand (with the complicity of 208 astrocytes 30 ) approximately 3-fold lower energies than their terminals 1,31 , albeit these estimations have been 209 calculated on the amount of mitochondrial ATP consumption. Accordingly, although authors reported that 210 oligodendrocytes, abundantly populating white matter 29 , provide metabolic support to neurons through 211 monocarboxylate transporters 32,33 , our results showed low vascular density in the white matter regions mostly 212 characterized by the presence of large vessels in terms of diameter. 213 The olfactory bulbs are important districts for the olfactory information processing, crucial for rodent behavior 214 and survival. The olfactory bulbs are the second neuron densest brain structure 29 (~250000 neurons per mm 3 ). 215 In olfactory bulbs tortuosity was very high (the second highest) indicating that vessels are highly twisted and 216 curved. Vessels in this region were characterized by the highest branch point density, low segment length and 217 large diameter. 218 The cerebellum is an evolutionary ancient three-layered brain section distinguished by the highest neuronal 219 density of about 830000 neurons per mm 3 . As expected, our estimations showed that cerebellar vessels were 220 the smallest in terms of diameter and one of the shortest according to segment length. Nonetheless, cerebellar 221 vessels had the highest distance to nearest vessel among regions, an unexpected result which suggested that 222 the ratio of vascular endothelial cells to the number of neurons is relatively low in comparison to other 223 regions 28 . 224 The striatum is a subcortical region responsible for motor functions with a relatively low neuronal density 28,29 225 (~64000 neurons per mm 3 ), and is generally divided in its ventral and dorsal (caudate putamen) parts. Ventral 226 striatum (~66000 neurons per mm 3 ) had the longest vascular segments and relatively low diameters, 227 compatible with a moderate level of energy demand. 228 At last, the neocortex, which is responsible of most of high-level brain functions and is characterized by a 229 moderate neuronal density (approximately 83000 neurons per mm 3 ), did not stand out in any vascular feature. 230 This result seems to be in accordance with recent 18FDG-PET measurements of the homogeneous metabolism 231 of the mouse cerebral cortex among other hindbrain and forebrain structures 34 . 232 As a summary, we presented a methodological framework for comprehensive and precise reconstruction of 233 the entire microvasculature of the mouse brain at the unprecedented pixel size of 0.65 µm. Local synchrotron-234 based X-ray phase-contrast tomography combined with an attainable computational pipeline resulted in an 235 effective methodology to investigate geometrical, morphological and topological features of vascular systems 236 of ex vivo organs at their finest structure. 237 The approach proposed in this paper has a few limitations. First, the whole brain sample analyzed in this work 238 was imaged in approximately 1200 tiles that required image acquisition session lasting more than 57 hours. 239 Each tile was therefore imaged in approximately three minutes. The sample must be steadily mounted and 240 stable such that during each three-minute interval it moves less than one pixel (0.65 µm), or otherwise the 241 tomographic reconstructions of the individual tiles may contain artifacts. Although this requirement could still 242 rise problems in various experimental setups, it is much easier to achieve than similar stability over the whole 243 imaging session. Note that deformations between neighboring tiles are acceptable in the stitching method used 244 here 15 . Second, the results shown here are based on a single animal and, although related literature does not 245 indicate important variations in the brain microvascular architecture, conclusions of this work could be slightly 246 different in a larger animal sample. Third, it seems to be hard to perfuse all vessels adequately with the 247 proposed contrast agent (Indian Ink), and therefore the non-perfused vessel branches are missing from the 248 analysis. This leads to biased results especially for the largest vessels, but according to results shown in 249 Section 2 and Supplementary Table 1 the smaller vessels seem to be unaffected. 250 In the future, the methods proposed in this work can be used, e.g., to construct an atlas of microvessel geometry 251 in mouse brain, both in healthy and pathological conditions, or to study blood flow in more detail using image-252 based flow simulations either in direct image-based modality 35 or using the generated vascular graph 36 . In 253 conclusion, we demonstrated that it is possible to make high-quality tomographic images of very large, fragile, 254 moisture-and radiation sensitive samples, and analyze their structure with image-based measurements. We 255 believe that the methodology introduced here generalizes well to many kinds of biological and engineered 256 samples, and is particularly useful in cases where optical clearing required in many other imaging modalities 257 is not possible or desirable. 258

Mouse sample preparation 260
The experimental procedure was approved by the local veterinary authority of Canton Zürich, Switzerland 261 (license number ZH184/2015). 262 After the loss of any reflex, before the death by the barbiturate overdose, the animals were prepared for the 263 perfusion by the opening of the sternal plate and the thoracic cage. The beating heart was then gently clamped 264 with flat tweezers and an 8-gauge metal needle with the smoothed tip was inserted into the left heart ventricle. 265 The intracardiac perfusion was performed in a few stages: first with a Ca+/Mg+-free phosphate buffered saline 266 (PBS) (100 ml, 37 °C), followed by a 4 % paraformaldehyde and Karnovsky's fixative (100 ml, 37 °C), and 267 then perfused with Indian Ink (50 ml, 40 °C), and finally with Karnovsky's fixative (20 ml, 4 °C). A clear sign 268 of the complete perfusion was the generalized blackening of all the mucosae, of nude skin surfaces (such as 269 the snout, the paws) the thoracic viscera and the eyes. form an image mosaic, where overlap between neighboring images is 30% of their diameter in directions 278 perpendicular to the rotation axis, and approximately 10% in the direction parallel to the rotation axis. 279 Each tomogram was reconstructed from 1001 X-ray projection images with the GridRec algorithm 37 . Paganin 280 phase retrieval method 38 was used before reconstruction. The projection images were acquired with 20 keV 281 monochromatic X-ray beam, 0.65 μm pixel size, and 50 ms exposure time. The sample to detector distance 282 was set to 100 mm. The total acquisition time was approximately 57 h. 283 284

Stitching 285
We stitched the individual tomograms into one large volume image using a non-rigid stitching algorithm 15 . 286 There, the locations and the orientations of the individual tomograms are globally optimized such that 287 disagreements between them in the overlapping regions are minimized. Furthermore, the overlapping regions 288 are deformed such that any remaining disagreements are eliminated. This processing ensures that the 289 microvessels are continuous across boundaries of individual tomograms even in cases where the sample has 290 deformed during image acquisition. The size of the stitched volume image was approximately 13600 × 13600 291 × 30000 pixels (11 TB with 16-bit pixels). 292 293

Registration to the Allen atlas 294
The stitched volume image was downsampled to similar size than the annotated Allen adult brain atlas at its 295 full resolution (version CCF 2017) 20 . The Atlas was then registered with the volume image, initially by an 296 affine transformation, and further refined with a non-rigid B-spline transformation. 297 298

Image segmentation 299
The stitched volume image was denoised using bilateral filtering 39,40 (spatial σ = 1.3 μm, radiometric σ = 7.6% 300 of full dynamic range), followed by high-pass filtering to remove large-scale intensity variations (spatial σ = 301 6.5 µm). The filtered image was segmented with a region-growing approach. To that end, initially all pixels 302 whose value were above a threshold were classified as vessels. The vessel regions were grown until all 303 bordering pixels had a value below a second threshold. All other pixels were classified as background. The 304 threshold values were selected such that the segmentation visually corresponded to the vessels. 305 Possible gaps in the segmented vessels were eliminated by applying a morphological closing filter (radius = 306 3.25 μm). In addition to blood vessels, the segmentation process identified choroid plexuses, some small and 307 separate non-vessel regions, and many planar structures (at the surfaces of the brain, mostly caused by 308 remaining phase contrast artifacts) as vessels. In order to eliminate the small non-vessel regions, all foreground 309 objects less than 685 μm 3 (equivalent to 2500 pixels) in volume were discarded. As the blood vessels form a 310 continuous network, this process does not have any effect on them. The choroid plexuses in the ventricular 311 cavities were removed by masking the original segmentation with a mask where choroid plexuses were not 312 visible. The mask was generated using morphological opening and closing operations. 313 Planar structures were eliminated by calculating a surface skeleton 41 of the foreground, and eliminating all 314 surfaces consisting of more than 5000 pixels. The surface skeleton was then refined into a line skeleton 41 where 315 each blood vessel is turned into a single pixel thick line located in the middle of the vessel. 316

Image analysis 318
The line skeleton was traced in order to produce a graph representation of the microvessel network. In the 319 graph, vessel bifurcation points are represented as vertices and vessels as edges. The bifurcation points were 320 found and the center lines between them, representing individual branches of the blood vessel network, were 321 traced in order to produce a graph representation of the microvessel network. For each bifurcation point, the 322 corresponding anatomical region was recorded based on the annotated volume registered with the image (see 323 Section "Registration to the Allen atlas"). For each microvessel branch, the length L of the branch 42 , distance 324 D between its end points, and the cross-sectional area A of the vessel was recorded. The cross-sectional area 325 was measured by taking two-dimensional cross-sectional slices of the vessel and measuring its area from 326 those 43 . The effective diameter of the vessel was then determined as = 2� ⁄ , tortuosity as / , and 327 slenderness as / . Branches shorter than 9.75 μm (equivalent to 15 pixels) and not connected to multiple 328 other branches in both ends did not correspond to vessels and were pruned. Finally, the distance between the 329 microvessels was quantified by calculating a distance map 44 where each non-vessel pixel is associated the 330 distance to the nearest blood vessel. 331 332

Uncertainty analysis 333
The uncertainty limits were estimated using a Monte Carlo method, where the image analysis process is 334 repeated several times with perturbed input parameters and the uncertainty limits are calculated from the 335 distribution of the results. In order to speed up processing, uncertainty analysis was done on 58 blocks of the 336 original volume image, 1500 3 pixels each, selected randomly from all anatomical regions in the brain. 337 The values of the input parameters were drawn randomly from normal distributions with means given by the 338 values used for analyzing the whole volume image, and standard deviations of 10% of the mean, except for 339 radiometric σ and threshold values where 5% and 2.5% were used, respectively.     Correlations between average neuronal (first and second column) and glial (third and fourth column) densities 516 and various measured quantities. Density data is from the EPFL mouse brain atlas 28