Dirigee Par

Nicolas Renier
email: nicolas.renier@icm-institute.org

Fallani Fabrizio De Vico

Devant

Serge Charpak

Sylvie Lorthois

Martine Cohen
email: martine.cohen-salmon@college-de-france.frspécialisteen

Salmon : Rapporteur
email: rapporteurlorthois@imft.frspécialisteen

Anatole Chessel
email: examinateuranatole.chessel@polytechnique.eduspécialisteenbio

Christoph Kirst
email: christoph.kirst@ucsf.edu

Sophie Skriabine

Alba Vieites-Prado

Nicolas Michalski

Marc Tessier-Lavigne

Thomas Topilko

Paul Bertin

Gaspard Gerschenfeld

Florine Verny

Piotr Topilko

Fabrizio De

Vico Fallani

Mapping of the mouse cerebral vascular network in adults and during development par Sophie Skriabine These de neurosciences

Binarization pipeline strategies, in particular works on adaptive thresholding and Hessian filters, DNN for vessel filling, and the development of new graph algorithms for the modeling and characterization of the vascular network.

Acknowledgement

I would like to thank Serge Charpak, Sylvie Lorthois, Martine Cohen Salmon, Anatole Chessel, Pablo Blinder, and Andy Shih for accepting to be part of my jury committee. This work would not have existed without the support and guidance of my supervisors, Nicolas Renier, and Fabrizio de Vico Fallani. I would especially like to thank them for giving me the opportunity to receive transdisciplinary supervision, be part of two amazing labs, and shaping me into the scientist I am today. I spent 5 fantastic years at the ICM. I especially thank Nicolas for welcoming me into his lab as an engineer before starting the PhD, for teaching me all about the biology of neuroscience, and for his scientific rigor, which allowed me to pursue a PhD under the best conditions. I am deeply grateful to all my lab mates from the RENIERLAB, particularly Elisa and Domitille, for their good spirits, support, and collaboration. Special thanks to Elisa, without whom this project would not exist, and to Domitille for teaching me so many new French expressions. I also thank Charly for his patience and valuable insights in coding and electronics, and Alba for her wise advice.

I also extend my gratitude to my lab mates from ARAMIS, especially Tristan Venot, Charley Presigny, and Vito Dicchio, for the running and climbing sessions, moral support, and weekly after-work gatherings.

I would also like to thank all our collaborators for their efficiency and dedication: Nicolas Michalski for the otoferlin mice and work on the auditory system, Alexandre Dubrac for the spatial transcriptomics, and Sophie Pezet for the FUS dataset.

I also thank Patricia Gaspar and Etienne Doumazane for proofreading this manuscript.

I give special acknowledgment to my parents who instilled in me a passion for science, especially my mother for her support, and my father, who would have been proud.

Finally, I thank Martin and all my friends for always being there.

Resume en Anglais

Cerebral vascularization is a vast network aimed at adequately distributing the nutrients and oxygen necessary for energy. Its functioning is closely linked to the activity of the neurons that populate it. This delicate balance involves multiple mechanisms, including neurovascular coupling interactions. It is assumed that neurovascular imbalances are involved in most brain and vascular diseases, such as neurodegenerative diseases. The precise mechanisms that control neurovascular plasticity remain largely unknown.

However, studying the role that the vasculature plays in all cognitive systems is a real technical challenge because very few resources are available. Therefore, we have developed a pipeline to image the vasculature of whole organs in a reasonable time frame in order to characterize the state and evolution of the vascular system under different pathological conditions [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF]. We also propose to monitor the evolution of the vascular network during development and provide resources in the form of atlases and computational tools to unlock research on brain development in mice.

This leads to the complete characterization of the vascular system: strong disparities were observed in both densities and topology. Particularly in the cortical areas, these differences are sufficient to characterize the function (integrative or somatosensory) of the different cortical regions. The observation of vascular plastic remodeling after stroke, as well as vascular degradation in congenitally deaf mice, provides strong evidence for the correlation between neuronal activity, stress, and vascular topology.

To f u r t h e r i n v e s t i g a t e v a s c u l a r h e t e r o g e n e i t i e s , w e e x a m i n e d t h e i r e v o l u t i o n d u r i n g d e v e l o p m e n t : Vascular development occurs in waves whose length and intensity depend on the brain region. The somatosensory systems (cortical and thalamic regions) seem to develop together, and it can be hypothesized that there is a link between the length of the developmental window and the orientation/density heterogeneities.

Early disruption of the system (e.g., by sensory deprivation) leads to a change in vascular topology in the adult, suggesting that decreased neuronal activity during development leads to decreased vascular density in the affected regions. In adulthood, such a loss of vessels has not been observed after a decrease in neuronal activity. However, a decrease in neuronal activity shows an impact during long-term aging in mice by accelerating vascular degradation, not only in somatosensory areas but also in thalamic areas related to memory and cognition. These areas are known to be involved in age-related neurological disorders in humans.

Resume en Francais

La vascularisation cérébrale est un vaste réseau visant à distribuer de manière adéquate les nutriments et l'oxygène nécessaires aux tissus cérébraux. Son fonctionnement est étroitement lié à l'activité des neurones qui le peuplent. Cet équilibre délicat fait intervenir de multiples mécanismes, notamment les interactions de couplage neurovasculaire. On suppose que des déséquilibres neurovasculaires sont impliqués dans la plupart des maladies cérébrovasculaires, telles que les maladies neurodégénératives. Les mécanismes précis qui contrôlent la plasticité neurovasculaire restent largement méconnus.

Cependant, étudier le rôle que joue la vascularisation dans l'ensemble des systèmes cognitifs représente un véritable défi technique car très peu de ressources sont disponibles. Nous avons donc développé un pipeline permettant d'imager la vascularisation de l'ensemble des organes dans un laps de temps raisonnable afin de caractériser l'état et l'évolution du système vasculaire dans différentes conditions pathologiques [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF]. Nous proposons également de suivre l'évolution du réseau vasculaire pendant le développement. Nous fournissons des ressources sous la forme d'atlas et d'outils informatiques pour debloquer la recherche sur le développement cérébral chez les souris.

Cela nous conduit à la caractérisation complète du système vasculaire: d'importantes disparités ont été observées tant en termes de densités que de topologie. L'observation d'une restructuration plastique vasculaire après un accident vasculaire cérébral, ainsi que la dégradation vasculaire chez les souris atteintes de surdité congénitale, apportent de solides preuves de la corrélation entre l'activité neuronale, le stress et la topologie vasculaire.

Pour approfondir l'étude des hétérogénéités vasculaires, nous avons examiné leur évolution au cours du développement : le développement vasculaire se produit par vagues dont la longueur et l'intensité dépendent de la région cérébrale. Les systèmes somatosensoriels (régions corticales et thalamiques) semblent se développer conjointement, et l'on peut émettre l'hypothèse d'un lien entre la durée de la période de développement et les hétérogénéités d'orientation/densité. Une perturbation précoce du système (par exemple, par la privation sensorielle) entraîne une modification de la topologie vasculaire chez l'adulte, suggérant que la diminution de l'activité neuronale pendant le développement conduit à une diminution de la densité vasculaire dans les régions concernées. La diminution de l'activité neuronale impacte également le vieillissement, en accélérant la dégradation vasculaire, non seulement dans le cortex, mais aussi dans les zones thalamiques liées à la mémoire et à la cognition, et reconnues pour être impliquées chez l'homme dans les troubles neurologiques liés à l'âge.

Introduction

"Throughout the body the animal arteries are mingled with veins and veins with arteries, and both veins and arteries are mingled with nerves and the nerves with these... And of course the usefulness of such a complete interweaving is very evident, if, that is to say, it is a useful thing for all parts of the animal to be nourished." --Galen, 2nd century A.D.

The vasculature has long been seen as a nourishing tree spanning across the body. However, until middle age, the relationship between veins and arteries was not completely understood, as veins were believed to carry blood and arteries carried spirit and energy through the body. Later anatomical studies persisted in presenting venous and arterial systems as completely separate, until perfusion experiments performed on vivisected animals showed that they belong to the same circulatory system (Harvey, 1660).

In the late 1800s, Angelo Mosso, Charles Roy, and Charles Sherrington [START_REF] Roy | On the regulation of the blood supply of the brain[END_REF] hinted at the possibility that brain activity increases cerebral blood flow. In the 1950s-60s, Seymour Kety and Lou Sokoloff developed autoradiographic methods to image regional cerebral blood flow during neural activation, jump-starting a novel field of research to understand neurovascular interactions. Indeed, due to their close cellular interactions within the neurovascular unit, the vasculature, astrocytes, and neurons strongly regulate each other's functions.

Without a blood supply, neurons start dying within 5 minutes. It has been demonstrated that the vasculature plays a major role in most, if not all, brain and systemic disorders, ranging from dementia to diabetes.

Even respiratory tract infections, such as those caused by SARS-CoV, unexpectedly, may also impact the cerebral vasculature in mice, thus increasing the risk of brain-related disorders.

However, studying neurovascular interactions is a daunting task, as both neuronal and vascular networks have a complex and very dense three-dimensional organization that spans several scales. Vessels vary in size by several orders of magnitude between capillaries and large veins, and neuronal connectivity can affect both local circuits and distant regions. Moreover, both networks are plastic structures that can be reorganized at small-and long-time scales. Obtaining reconstructions of the organization of both networks represents an imaging challenge because one needs to image centimeter-sized regions with micrometer resolution.

Additionally, owing to the variability, both in time and across individuals, of the topology of these networks, scans of several brains are needed to derive robust information on their organizational principles. Another challenge is the identification of the types of vessels in the imaging datasets (arteries, capillaries, veins).

Therefore, imaging data capturing the organization of the vascular network is of a very large size, which calls for the development of novel, tailored image processing strategies. Another challenge to address is that, even though the vascular system can be formally represented as a network from a mathematical point of view, very few algorithms can scale to the number of nodes and edges needed to embed the full brain vasculature.

Moreover, the spatial heterogeneity of the cerebral vascular network may require different analysis tools depending on the brain regions considered. To summarize, understanding the structure of the cerebral vascular network requires:

• Reproducible and parallelized brain preparations (i.e., avoiding the need for complex perfusions)

• Acquisition systems compatible with large intact volumes (i.e., above 1 cm3)

• Adapted spatial resolution to capture brain capillaries (i.e., smaller than 3 µm/pixel)

• Fast data acquisition to obtain many individual scans

• Annotation of the vessels: recognizing veins, capillaries, and arteries

• Fast algorithms that scale to the image sizes (i.e., capable of processing over 100 GigaPixels)

• Fast algorithms that can scale to the size of the graphs for analysis (over 100 million nodes)

The major roadblock up to this day to obtain complete reconstructions of the vascular network has been limitations in both the availability of high-performance algorithms and computers, as well as sample preparation techniques that enable streamlined 3D acquisitions. The few datasets available today were obtained after very intensive efforts and constitute the foundations of a novel field modeling blood flow in silico. However, since those initial efforts, there hasn't been much progress done until the recent advances in tissue clearing and light sheet microscopy.

In this thesis, I aimed to develop tools for the reconstruction and analysis of oriented graphs (arteries->capillaries->veins) of the cerebral blood vessels based on light sheet microscopy and tissue clearing.

Combined with mouse genetic tools developed by our collaborator Nicolas Michalski at the Pasteur Institute/Audition Institute to manipulate sensory experience and neuronal activity in the auditory system, I strived to generate and analyze the vascular graphs of mouse brains obtained from lines with congenital or progressive deafness, as well as deaf mice having undergone hearing restoration. To further develop the tools and to figure out how those networks are impacted by neuronal plasticity, I also generated the resources to study the postnatal developing brain with the generation of an MRI/Lightsheet postnatal developmental anatomical atlas.

My project aims to achieve the following challenges:

• Map the topology of vascular networks through their interindividual variations

• Define the interactions between neuronal activity and vascular plasticity during development

• Develop novel algorithms and methods aimed at better understanding complex biological graphs containing millions of nodes

• Understand the interactions between the heterogeneous metabolic demand of the brain and the vascular topology

In the first part, I will give a brief overview of the basic biology of neurovascular coupling. In the second part, I will present and discuss state-of-the-art techniques, both in terms of microscopy and computational pipelines, for the generation of 3D vascular datasets. This part aims to pinpoint the pros and cons of the different existing pipelines and justify how the method I contributed to developing improves the current state of the art and outperforms the existing applications. In the third and last part, I will introduce the concept of graphs, which is the format used to represent and extract the vascular system. I will present the classical methods for traditional graph analysis and explain how they are limited in the case of the vascular architecture. Then, I will describe the vascular topology as a graph and present the alternative strategies to make the most out of the graph format by extracting different information, and how it makes room for novel algorithms able to tackle large graph data.

Then I will present in greater details the work that has been done during the thesis, and finally I will discuss the limitations and prospecting improvements of different aspects of the project.

I. The vasculature, a supply highway

The brain vasculature aims to distribute blood, oxygen, and nutrients to adequately fuel the neural cells, especially the neurons. The scaling of neurons across mammalian species varies linearly: as the number of neurons increases, neuronal density decreases. However, due to the scarcity of available data, little is known about how the vasculature itself scales in comparison to neurons and how it affects its structure (Fig 2).

I.1 General structure of the vasculature

The brain is a crucial, yet fragile and dynamic organ that undergoes continuous and subtle activity modifications over time (Iadecola, 2017). It relies on a precise and adaptable energy supply for its proper functioning. To achieve this, it is perfused by a dense and intricate network of arteries, capillaries, and veins.

The complex interaction between the neuronal and vascular networks enables subtle changes in neuronal activity to influence local blood flow in the mature brain, as well as vascular development and remodeling in the developing brain [START_REF] Carmignoto | The contribution of astrocyte signalling to neurovascular coupling[END_REF].

Although the brain accounts for roughly 20% of the body's total energy consumption, it has limited energy storage capacity. Therefore, it heavily depends on a continuous and stable blood supply from the vascular network to receive energy, oxygen, and nutrients. Blood is supplied to the brain through the carotid arteries, which originate from the aorta, and the vertebral arteries, which connect to the subclavian artery. The carotids connect to the circle of Willis, which distributes blood to the major cerebral arteries of the forebrain, while the vertebral arteries feed into the basilar artery of the brainstem. Notably, the surface arteries of the brain have anastomosis, providing alternative routes to compensate for potential strokes or arterial blockages. Pial arteries run along the surface of the parenchyma, giving rise to penetrating arteries, arterioles, and eventually capillaries. Blood is then drained into venules and finally into veins [START_REF] Carmignoto | The contribution of astrocyte signalling to neurovascular coupling[END_REF].

Pial veins collect blood from the superficial parenchymal regions, while a large venous central sinus collects most of the blood from the deep cerebral regions. The gross anatomical structure of the mouse cerebral vasculature has been accurately reconstructed using the fMOST (fluorescence micro-optical sectioning tomography) technique.

In the human brain, which has a volume of approximately 1200 cm³, the total vascular network contains up to 640 km of capillaries (Pandey et al., 2016). In mice, the capillary network reaches 384m for a cranial volume of approximately 1.5 cm³ [START_REF] Ji | Brain microvasculature has a common topology with local differences in geometry that match metabolic load[END_REF]. Rather than a tree-like graph where capillary branches divide sequentially from the arteries, the vascular network resembles a three-dimensional web with highly redundant paths. Its redundant and interconnected structure allows for efficient nutrient distribution and provides a certain degree of protection against strokes by allowing blood flow to be reversed in case of occlusion [START_REF] Schaffer | Two-photon imaging of cortical surface microvessels reveals a robust redistribution in blood flow after vascular occlusion[END_REF].

I.2 Cellular structure of the vasculature

Blood vessels consist of endothelial cells, which line the inner part of all vessels and originate from the bone marrow in adult animals [START_REF] Zhao | Establishment and Dysfunction of the Blood-Brain Barrier[END_REF]. In the brain, these endothelial cells undergo a specific specialization program that downregulates the general trans-cellular machinery responsible for molecule transport. Instead, they develop selective transport mechanisms. Additionally, brain endothelial cells express high levels of tight junction proteins, such as Claudin 5, which restrict the movement of solutes around the cells. Consequently, these cells form the blood-brain barrier (BBB), the primary component that highly restricts the permeability of the brain vasculature [START_REF] Daneman | The blood-brain barrier[END_REF].

The limited permeability of the BBB prevents rapid changes in ionic or metabolic conditions by restricting the movement of substances from the systemic circulation to the brain. This buffering effect helps maintain a stable environment for neuronal function. Moreover, the BBB protects the brain from exposure to substances that are harmless to peripheral organs but toxic to brain neurons, such as glutamate, nitric oxide, or botulinum toxin (i.e., botox). The permeability of the BBB is influenced not only by neurons but also by the extracellular matrix and non-neuronal cells, including astrocytes, pericytes, and vascular endothelial cells [START_REF] Sweeney | The role of brain vasculature in neurodegenerative disorders[END_REF].

Astrocytes are star-shaped glial cells contacting both neurons and blood vessels. One of their many roles is to mediate the passage of nutrients from the bloodstream to the surrounding neurons, regulate the microchemical environment, and modulate neuronal activity. Most, if not all, astrocytes have at least one endfoot unsheathing a blood vessel. [START_REF] Armulik | Endothelial/pericyte interactions[END_REF] Pericytes are cells that are wrapped around vessels, particularly the capillaries. They are in contact with several endothelial cells with whom they communicate via specialized peg-socket connections. Therefore, they are thought to play a role in the communication and coordination of endothelial cells within their environment. They are also known to be involved in vasculogenesis, angiogenesis, and can control the opening of the blood-brain barrier [START_REF] Winkler | Central nervous system pericytes in health and disease[END_REF].

Endothelial cells, astrocytes, pericytes, and neurons, along with their extracellular matrix, function as a neurovascular unit and regulate BBB permeability and maintain the integrity of the central nervous system. However, under certain conditions such as inflammation, traumatic brain injury, or ischemic stroke, the BBB can be compromised, allowing for the passage of larger hydrophilic substances [START_REF] Lacoste | Control of cerebrovascular patterning by neural activity during postnatal development[END_REF].

Cerebral blood vessels larger than capillaries are lined with mural cells. A diversity of pericytes and smooth muscle cells wrap around the endothelial cells, with cellular morphologies and molecular signatures that gradually change in the artery -> capillary -> vein axis, as a continuum. Pericytes control blood flow in capillaries and are important in many biological processes, and their role in diseases has been more recently recognized. Arteries and arterioles are covered with an extra layer of smooth muscle cells organized as rings, allowing tight adjustments of the intraluminal vessel diameter that can also regulate the blood flow. [START_REF] Hillman | Coupling mechanism and significance of the BOLD signal: A status report[END_REF] The regulation of cerebral blood flow (CBF) itself is still a matter of intense research. Over the past years, various hypotheses concerning the possible roles of the diverse abovementioned cell types (pericytes, astrocytes, smooth muscle cells, endothelial cells) have been proposed and debated. While it seems clear that the regulation of the CBF can occur at different levels, neuronal activity is thought to back-propagate from the site of increased demand all the way to the pial surface and hence act on several cell types. Pericytes, acting on capillaries, may have the strongest effect on CBF, as they are positioned at the site of the highest resistance in the vascular network (Fig 3). Nevertheless, we still do not have an integrated mechanism describing the coupling between neuronal activity and CBF across the different cell types (see table 1 for a summary of the known actors).

Table 1: from E.M.C Hillman 2014 : the table below summarizes the known signaling pathways bridging neuronal activity with vasomotricity, with a summary diagram presented below. Most neuronal and glial cell types interface with neuronal activity via diverse molecular mechanisms to modulate the diameter of capillaries and arterioles through the contraction of pericytes and smooth muscle cells.

Cell

I.3 The cerebral vasculature in diseased states

Because of its tight integration with neural cells, the vasculature almost always reacts to some extent to brain diseases and is, for instance, the focus of intense research in the context of neurodegenerative disorders.

Sometimes, the vasculature is an aggravating factor, but not the cause of an ailment, while other times it is the direct cause of brain function impairment.

A

For instance, traumatic brain injuries (TBI) can illustrate how secondary alterations to the brain vasculature may be the cause leading to secondary neurodegeneration long after the injury, even in the case of mild initial injuries. TBIs are an extremely common and widespread brain affliction. In the U.S, 1.7 million TBIs occur each year, ranging from mild to severe. Most of the time, they are the result of repeated shocks to the head (sports), car accidents, or falls, and have short-term consequences. However, even mild TBIs have been recently associated with cognitive deficits (memory or reasoning loss) and abnormal social behavior (anxiety or poor social skills). Additionally, they increase the risk of later developing dementia. Chronic damages such as vasospasm, BBB weakening, and chronic inflammation are among the pathologies associated with traumatic brain injuries (TBI). Studies in rodents show that most vascular damage can be healed within 2 to 3 weeks after the TBI. However, even if the vascular density is set back to baseline, sub-optimal cerebral blood flow restoration is observed, revealing widespread hypoperfusion regardless of the structural recovery of the vascular network (Fig 4). Furthermore, the newly formed vessels after severe strokes presented morphological abnormalities such as increased tortuosity and impaired BBB, but also vessel wall aberrations and reduced vessel radii. Both capillaries and larger vessels seem to be affected, and these deficits lasted longer than a few months after the injury. This suggests that even though vessel density can be restored, TBI can lead to delayed hypoperfusion and increase local damages due to poor vascular repair.

Peripheral metabolic syndromes can also directly affect the cerebral vasculature. For instance, the vascular factor in diabetes has long been described. It is a widely accepted fact that patients affected by diabetes are more likely to develop nephropathy, atherosclerosis, and retinopathy since hyperglycemia is associated with hypertension, dyslipidemia, and atherothrombotic complications. [START_REF] Rask-Madsen | Vascular complications of diabetes: mechanisms of injury and protective factors[END_REF]. However, in addition to hyperglycemia alone, the side effect of insulin resistance seems to play a role in the development of vascular abnormalities. Most of the pathologies linked with diabetes are due to inflammation processes and apoptotic reactions.

Vascular dysfunction is increasingly recognized as a key component of neurodegenerative disease (ND) as it increases the risks of dementia and early cognitive decline. While genetic, lifestyle, and environmental factors have been demonstrated to be correlated with AD, they also influence neurovascular health. In addition to its An outstanding question is what drives the burst of neovascularization during the postnatal first 2 weeks.

Among all the possible candidates, myelination and synaptogenesis occurring simultaneously with the angiogenesis postnatal burst constitute likely actors in this phenomenon. Axonal branching and dendritogenesis could transiently increase hypoxia by shaping neuronal activity and its associated metabolic cost, thus driving angiogenesis.

Microglia invasion also happens early during embryonic development [START_REF] Thion | On place and time: microglia in embryonic and perinatal brain development[END_REF]. From E16, microglia start colonizing neocortical areas, modulating many developmental processes such as synaptic plasticity. Since they are also closely associated with blood vessels [START_REF] Goldmann | Origin, fate and dynamics of macrophages at central nervous system interfaces[END_REF], they may directly or indirectly play a role (by modulating other cells' growth and oxygen needs) in angiogenesis during development [START_REF] Fantin | Tissue macrophages act as cellular chaperones for vascular anastomosis downstream of VEGF-mediated endothelial tip cell induction[END_REF][START_REF] Rymo | A two-way communication between microglial cells and angiogenic sprouts regulates angiogenesis in aortic ring cultures[END_REF].

Neuronal activity is also a factor controlling postnatal angiogenesis. Lacoste et al. (2014) showed that somatosensory deprivation in one-week-old pups leads to a decrease in endothelial cell proliferation and, therefore, a sparser and less branched vascular network in the somatosensory cortical areas. Conversely, they also showed that increased sensory input had the opposite effect. Oligodendrocytes, responsible for myelination, also mature concomitantly with blood vessels and closely interact for their migration and for the stimulation of angiogenesis (Tsai et al., 2016). Polydendrocytes, NG2+ cells, are also necessary for the sprouting of new vessels between E12 and E18 and likely also in postnatal stages [START_REF] Minocha | NG2 glia are required for vessel network formation during embryonic development[END_REF].

In conclusion, vessel formation seems to result from a complex and intricate mix of interactions among various actors. While most of the molecular actors are known, the coordinated maturation of the system in the postnatal stages and the integrated role of each cellular player are still being refined. Access to an efficient methodology for the complete 3D reconstruction of the network would greatly assist in determining the precise contribution of each cell type to the mature vascular structure.

II. Generation of cerebral vascular datasets

As mentioned in the introduction, methods to acquire the 3D structure of the cerebral vascular network are necessary to better understand its organization, development, and plasticity. It is challenging to design a method that can simultaneously capture the complexity of the vascular network while still being as automated and fast as possible. The brain is an incredibly complex organ, composed of millions of delicate and fragile vessels with dimensions and length ranges spanning 1 or 2 orders of magnitude. Additionally, there can be regional heterogeneities in terms of density, orientation, and tortuosity. Therefore, high-resolution 3D acquisition at a large scale is required. To provide an extreme example, classical histological techniques combined with fluorescent microscopy work at subcellular resolution but only for micrometer-thick slices. On the other hand, MRI allows fast 3D acquisition of whole-brain data but has limited spatial resolution. Another challenge lies in the processing of the acquired data. Imaging the entire organ at the capillary level results in hundreds of gigabytes of data, which requires significant computational power to be processed within a reasonable amount of time. Novel approaches based on machine learning networks, parallel computation, and the use of clusters have become key in bioimaging pipelines to address this challenge and have gained popularity.

II.1 Legacy methods for 3D vascular reconstructions

II.1.1 Resin casts

Since the invention of non-viscous polymer 60 years ago [START_REF] Murakami | Application of the Scanning Electron Microscope to the Study of the Fine Distribution of the Blood Vessels[END_REF] , vascular cast corrosion has been widely used to describe the normal and pathological vascular system of various organs. The technique involves injecting Mercox resin into the vascular system and then removing the surrounding tissues with KOH, hot water, and formic acid. The resulting 3D cast can be imaged using scanning electron microscopy [START_REF] Gaudio | A Tridimensional Study of Microcirculation in Skeletal Muscle[END_REF] . However, the assumption that the corrosion cast technique faithfully represents the entire capillary network is not completely correct. It has been shown that due to their delicate structure, some capillaries can collapse during the maceration phase, and some may not be filled with resin due to insufficient injection flow.

Additionally, Mercox resin is known to shrink during polymerization, ranging from 5% to 20% depending on the degree of dilution (Weiger et al., 1986). Moreover, the viscosity of the resin itself is highly sensitive to aging and external conditions, and it varies from one batch to another, making measurements hard to reproduce [START_REF] Hossler | Vascular Corrosion Casting: Review of Advantages and Limitations in the Application of Some Simple Quantitative Methods[END_REF] . The biggest drawback of this technique is that it produces a 2D image obtained from the resin surface via scanning electron microscopy, which has a limited field of view and a limited depth of imaging restricted to the sample surface.

To a d d r e s s t h i s i s s u e , m i c r o -computed tomography (micro CT) can be used to scan the resin cast. MicroCT utilizes tomographic reconstruction and non-invasive techniques such as 3D X-ray diffraction to produce 3D image acquisitions. This technique is widely used in biomedical imaging, similar to CT scanners. It enables relatively fast scanning of large samples and can achieve sub-micrometer resolution. For example, it takes approximately 15 hours to scan an entire mouse brain vasculature with a 2um voxel size, resulting in a 500GB dataset per sample (Quintana et al., 2019) . However, since this technique needs to be combined with the corrosion cast injection methods described above, it is prone to inconsistencies due to resin quality. The main drawback of this technique is that it is not possible to image the background brain tissue since it has been digested away, rendering the data unalignable to a reference brain.

II.1.2 Serial tomography

To o v e r c o m e t h i s l i m i t a t i o n , o t h e r t e c h n i q u e s a r e b a s e d o n s e r i a l t o m o g r a p h y. I n t h i s t e c h n i q u e , t i s s u e s a r e sliced using a cryostat. The slices are then mounted, stained, and imaged using microscopy. The images are then assembled. Since the method involves slicing, the axial resolution is limited by the thickness of the slice, introducing some anisotropy and distortions that prevent accurate 3D reconstruction. The technique is also labor-intensive, and the cutting process before imaging inevitably introduces deformations [START_REF] Amato | Whole brain imaging with serial two-photon tomography[END_REF]. I n e a c h iteration, a section is removed with a vibratome, and the fresh brain surface is imaged (Ewald et al., 2002).

However, even if some of the previously mentioned limitations are overcome, there are still significant drawbacks. Firstly, slicing introduces surface damage to the sample. Moreover, the slice thickness continues to limit the axial resolution. Originally, this technique was associated with whole-field imaging, so the intrinsic fluorescence of the underlying tissue adds noise to the imaged surface. The entire process is also timeconsuming. These issues have been addressed by using a laser to cut the surface layer in each iteration [START_REF] Winiarski | Broad ion beam serial section tomography[END_REF] . Another important improvement is the combination of two-photon microscopy with serial tomography [START_REF] Amato | Whole brain imaging with serial two-photon tomography[END_REF].

Two-photon microscopy is a high-resolution 3D imaging technique that relies on the simultaneous absorption of two photons, whose combined energy allows excitation of an electronic state. The fluorophore then relaxes to the first electronic state, releasing an electron. Due to the low probability of simultaneous electron absorption, which is highest near the focal point, two-photon microscopy has a relatively shallow Point Spread Function in the axial direction and sub-micron resolution in the imaging plane (Peter T.C. So et al., 2000) . In 2013, Blinder et al. used this technique in combination with fluorescent gel injections to image the vasculature of the mouse somatosensory cortex (Blinder et al., 2013) . They imaged a cube of 1900x1200x1450 um with a 1-micron resolution in the barrel field and segmented the vessels while manually labeling veins and arteries.

The data quality is excellent, and the overall process is mostly automated with simpler sample preparation.

However, it still has some drawbacks. First, it relies on the quality of dye injection, which can be affected by poor perfusions. More importantly, only a small volume can be imaged at a time. However, this approach has recently been scaled up to the whole brain.

A possible way to overcome these limitations is to work with histological staining, for instance, by imaging successive stained brain slices using two-photon microscopy [START_REF] Ragan | Serial two-photon tomography for automated ex vivo mouse brain imaging[END_REF] . The use of clearing methods combined with lightsheet microscopy would streamline this type of preparation, albeit with a loss in imaging resolution [START_REF] Vieites-Prado | Tissue clearing and 3D imaging in developmental biology[END_REF] .

II.2 Light sheet microscopy and tissue clearing

II.2.1 Light sheet microscopy

Since transparency has always been a highly valuable feature for observing organisms such as zebrafish and cells, techniques to optically clear tissues have been developed early on. At the beginning of the 20th century, Werner Spalteholz published a clearing protocol based on methyl salicylate and benzyl benzoate [START_REF] Spalteholz | A Method for the Clearing of Human and Animal Specimens[END_REF] . In 2007, the combination of clearing with light sheet microscopy paved the way for large sample and whole organ imaging and analysis [START_REF] Dodt | Ultramicroscopy: Three-dimensional visualization of neuronal networks in the whole mouse brain[END_REF] . Light sheet scans, being performed on intact samples, do not exhibit mechanical distortions and are well suited for acquiring large samples. However, the chemical tissue processing during optical clearing can cause some distortions.

II.2.2 Tissue clearing

Clearing is one of the crucial steps of whole brain imaging. It must be chosen carefully depending on the application, the type of organs, or the structures that will be targeted. Even though the different techniques vary in terms of chemicals, the overall idea is the same: the sample is first fixed to maintain its structure throughout the processing. Then, the sample is depleted of lipidic molecules to allow the clearing solution and antibodies to diffuse through the tissue. Then comes immunolabeling to stain specific structures, and lastly, the sample is flushed with an optical clearing solution to homogenize its refractive index [START_REF] Vieites-Prado | Tissue clearing and 3D imaging in developmental biology[END_REF] . Optical clearing is one of the crucial steps of whole brain imaging with light sheet microscopy. A clearing protocol must be chosen carefully depending on the application, the organ, or the structures that will be targeted. (Fig 10) Even though the different techniques vary in term of chemicals, the overall idea is the same:

II.3.1 The differents tissue clearing methodologies

the sample is first fixed to maintain its structure through the processing. Then, the sample is depleted of lipids to allow the refraction-index matching solution to homogenize the refraction index to the one of the proteins, making the sample transparent. Lipid clearing is also important in case of antibody labeling, as it is necessary for antibodies to diffuse through the tissue [START_REF] Vieites-Prado | Tissue clearing and 3D imaging in developmental biology[END_REF] .

Tissue clearing methods can be classified into three main families: 1) Aqueous based, 2) hydrogel-based, and

3) organic solvent-based methods.

Aqueous methods, such as CUBIC (Susaki et al., 2020) , tend to be simpler and non-toxic. They rely on hydration to wash away lipids as micelles. However, the protocol must be adjusted depending on the labeling and the type of sample and presents specific imaging challenges due to the sample's expanded size and lack of rigidity.

On the other hand, hydrogel-based methods such as CLARITY [START_REF] Treweek | Whole-body tissue stabilization and selective extractions via tissue-hydrogel hybrids for high-resolution intact circuit mapping and phenotyping[END_REF] , are designed to increase tissue stability and enable the use of stronger detergents. However, it is a very complex protocol with many variables to optimize depending on the desired application. The clearing method selected will depend on the size of the sample, the microscope used, and the nature of labeling (endogenous, tracers, or antibodies). For my applications of imaging the cerebral vasculature, iDISCO combined with a light sheet microscope of the Dodt design (Ultramicroscope) was one of the most convenient options.

Among the solvent-based techniques, many of them were designed to be associated with fluorescent proteins, with varying levels of success. However, a breakthrough was to realize that antibodies could diffuse deep in these types of heavily permeable samples. Therefore, the tissue can also be labeled intact using immunolabeling techniques instead, making the complex preservation of endogenous fluorescence unnecessary. Therefore, iDISCO has been developed to enable deep antibody diffusion in whole organ samples while keeping the processing simple enough to allow for parallel sample preparations with high throughput (Renier et al., 2014) . With standard PFA perfusion, it makes a cheap, reproducible, and reliable platform for the analysis of brain vascular networks. This method is:

• Not perfusion dependent

• No mechanical deformation due to slicing (although chemical distortions due to clearing)

• Micron resolution with the light sheet microscopy (1.6µm/pixel at 4X with 0.35NA)

• Compatible with additional histology

• Cheap to scale up: the cost for 1 brain is $20

• Highly transparent samples

• Light sheet compatible
The protocol includes the following steps: samples are fixed with PFA and harvested, and then washed in PBS.

They are then dehydrated with methanol, bleached, and then rehydrated. They are then immunolabeled using permeabilization and blocking solutions before being incubated in primary and then secondary or conjugated antibodies. Samples are then washed and cleared with DBE before being imaged with the light sheet microscope.

Volume imaging on light sheet microscopes generates extremely large volumes of data (200GB for one single channel of half a mouse brain). Those images must be processed, and the labeled vessels must be segmented efficiently. However, most computers do not have the computational power and RAM to treat such a dataset in a reasonable amount of time. Moreover, classical computer vision approaches do not necessarily scale to this type of large dataset. Therefore, novel pipelines must be carefully crafted and optimized in order to overcome this challenge .

II.3 Processing of light sheet generated images

Once the data is acquired, vessels must be reconstructed from the scans and analyzed. A faithful and accurate vessel reconstruction is essential to provide a means to reconstruct the vascular topology, with the goal being to estimate cortical blood flow in silico [START_REF] Blinder | Topological basis for the robust distribution of blood to rodent neocortex[END_REF] , (Schmid et al., 2017) , assess relationships between neuronal and vascular density, perform topological analysis (Blinder et al., 2013) , or witness vascular remodeling. However, an exact reconstruction is extremely challenging to obtain: From the sample preparation to the data acquisition and extraction, many errors can be made. For instance, if the sample harvesting is based on intracardiac injection, a trade-off must be found to ensure that even the smallest finest capillary is filled with the gel without breaking them. Then, the fixation and clearing pipeline introduce deformations and tissue shrinkage that must be minimized. In the case of serial two-photon tomography, the slicing of the sample itself can introduce artifacts. During imaging, artifacts can also occur. For instance, the signal-to-noise ratio should be as homogeneous as possible, but in the case of 3D antibody staining, it is known that diffusion through the tissue makes the signal in the center of the sample less good than on the surface. To correctly compute the flow resistance of a vessel, its radii need to be carefully measured, meaning acquisition resolution should be at a micrometer at least, or even less. Depending on the acquisition methods, reconstruction strategies might differ to address the challenges inherent to the techniques and correct specific artifacts.

However, all segmentation pipelines can be seen as modular toolboxes where each step can be modified or adapted to specific datasets. Typically, it can be divided into three large main modules that will be detailed underneath: pre-processing, segmentation, and post-processing (Fig 11).

II.3.1 Common Preprocessing steps

Pre-processing aims to make the data as clean as possible before the segmentation process. In most applications, it is impossible to rely solely on segmentation and post-processing operations. The main idea is to increase the distinguishability between foreground and background pixels to facilitate later segmentation, either by increasing the homogeneity within each class or the heterogeneity between classes (Arganda-

Carreras & Andrey, 2017).

Noise attenuation:

Noise can be due to tissue autofluorescence interfering with signals, unbound antibodies, or acquisition readout. It can impact the distribution of the foreground and background signal and lead to segmentation errors. The standard ways to address it are linear filters (Gaussian blur, mean filter), which also affect local contrast. Median filters, sigma filters, bilateral filters, or anisotropic diffusion are popular alternatives to avoid boundary blurring, but they come with the cost of having more parameters to tune.

Background correction:

Estimation of clipping values based on conservative considerations: the lowest voxel value in the background in the center of the brain and the largest on the periphery.

Contrast enhancement: Traditionally done using Laplace filters or deconvolution methods to remove blur.

Artifact correction: For instance, light sheet introduces horizontal striping artifact due to side illumination.

Stripes can be detected using Fourier analysis or by estimating the local intensity of the stripe according to the intensities of the voxel along the stripe artifact axis and an estimation of the local background intensity.

Deconvolution: Removes halos and blur introduced by the microscope optics or counteracts the effect of noise attenuation algorithms. Deconvolution procedures are computationally demanding and can be replaced by subtracting estimated high pixel halo values.

Object/geometrical specific enhancement: When the object of interest presents specific morphological features, these can be enhanced. Composite opening-closing operations or top-hat transform are classical ways of achieving this. For vessels and cell detection, 2nd derivative-based Hessian filters can be used to enhance their structure.

II.3.2 Segmentation steps

Before analysis, objects must be classified as background or foreground pixels. There are various methods to do so, and they can strongly depend on the dataset or data specifics. Among the most popular ways to achieve this are:

Binarization thresholding: It consists of comparing the value of a pixel to a threshold value. Any pixel above the threshold will be set to 1, and any pixel below will be set to 0. This results in a binary image discriminating the background from the foreground. However, especially in large 3D datasets, some regions might undergo unpredictable intensity variation, disabling the use of simple thresholding. Automated adaptive methods such as Otsu methods or maximum entropy enable compensation for region signal intensity variability.

II.3.3 Common post-processing steps

Even after all the previous steps, data might need refinement. For instance, for vessels stained with antibodies targeting endothelial cell proteins, the lumen can remain hollow. This may lead to topological aberrations when the user tries to extract the centerline of the vessel, for instance. Most corrections needed at this stage involve regularization of object shape (smoothing), filling holes (false negatives), and removing over-segmentation in the background (false positives). Dilation and erosion morphological filters are typically used in this case.

II.4 Deep learning approaches

II.4.1 Introduction to deep learning

Recent breakthroughs in imaging allow biologists to collect high-resolution data, resulting in hundreds of gigabytes of data. Therefore, computational power is now an essential tool. In the meantime, deep learning has become increasingly performant, outperforming traditional computer vision algorithms. Traditional approaches, such as the one described above, are handcrafted and therefore limited by the accuracy of the coder's assumptions. That is why data-based or statistical approaches might be more suited to solve this type of task. Convolutional neural networks (CNN) were particularly relevant in the field of image interpretation and analysis. However, despite its remarkable success, the technique comes with its share of challenges and limitations (Belthangady & Royer, 2019). In 2012, Krizhevsky et al. won the ImageNet Large-Scale Visual Recognition Challenge with an unprecedented margin, with the now-famous AlexNet. Since then, Convolutional Neural Networks (CNN) continued to break records and even surpassed human accuracy in 2018 [START_REF] Krizhevsky | Cnn实际训练的[END_REF] .

Deep learning allows computational models that are composed of multiple processing layers to learn representations of data with multiple levels of abstraction. It discovers intricate structures in large datasets by using the backpropagation algorithm to indicate how a machine should change its internal parameters that are used to compute the representation in each layer from the representation in the previous layer [START_REF] Lecun | Deep learning[END_REF] . It aims to transform an input into the output by passing it through a nonlinear parametrized and differentiable block module (layers). Their parameters are set during the training phase where each module is iteratively updated to improve the accuracy of the output to ground truth matching. It typically uses the backpropagation algorithm to do so [START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF] . Backpropagation uses the chain rule to efficiently compute derivatives and gradients. Each network aims to minimize a function, the loss, that computes the discrepancy between its output and the ground truth [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF] . The loss can be seen as a highdimensional landscape where gradient descent tries to find the lowest point. Convolutional neural networks are one of the most popular deep learning approaches when it comes to image analysis and reconstruction.

They can theoretically approximate any function if given enough parameters [START_REF] Masci | stacked convolutional AEs for hierarchical feature extraction[END_REF] .

II.4.2 Convolutional neural networks

A convolutional neural network is a network in which the connections between the nodes do not form a cycle.

It often contains up to 20 or 30 layers. The power of a convolutional neural network comes from a special kind of layer called the convolutional layer. Convolutional neural networks contain many convolutional layers stacked on top of each other, each one capable of recognizing more sophisticated shapes. With three or four convolutional layers, it is possible to recognize handwritten digits, and with 25 layers, it is possible to distinguish human faces.

The usage of convolutional layers in a convolutional neural network mirrors the structure of the human visual cortex, where a series of layers process an incoming image and identify progressively more complex features [START_REF] Masci | stacked convolutional AEs for hierarchical feature extraction[END_REF].

II.4.3 Common image segmentation architecture

Since input and output images in microscopy can be both 2D and 3D, architectures like the encoder-decoder are one of the simplest that can manage such an image-to-image transformation in multiple dimensions.

First, the input image is successively and repeatedly downsampled, convolved, and subjected to other normalization and nonlinear operations typical for CNNs. At each step, the image dimensions are reduced while the number of image channels-an additional dimension-is increased. This creates a representation bottleneck that is then reversed with successive and repeated upsampling and convolution steps. The bottleneck forces the network to encode the image into a small set of abstracted variables (often called latent variables) along the channel dimension, which is then decoded by the second half of the network. This ensures that only the most relevant features from the input image are retained, and those that are deemed unessential to the reconstruction are discarded. In the GAN (Generative Adversarial Networks) framework, two networks, a generator and a discriminator, are trained together. In microscopy applications, an improvement on GANs is often used: cGAN (conditional generative Adversarial Network) (Fig 12). The generator learns to transform an input image into an output image, while the discriminator learns to classify images as real (from the training set, label 1) or counterfeit (from the generator, label 0). The input image is provided to the discriminator. Therefore, cGANs can be understood as providing an implicitly learned loss function [START_REF] Goodfellow | Generative adversarial networks[END_REF] .

II.4.2 Limitations

II.5 3D Vascular reconstruction strategies

Given all the previously described tools there are many ways to create a full pipeline aimed to extract the vascular connectome at the whole brain scale with microcapillary precision. In the following paragraphs, I will describe a few of the most popular ones and compare their respective strength and weaknesses.

II.5.1 Serial 2p tomography -Ji et al., 2021

In this first study, no clearing was used. They perfused mice to replace the blood in the vessels with a fluorescently labeled gel. The brain is later dissected, fixed, and imaged with serial two-photon tomography at a voxel resolution of 0.3 * 0.3 * 1.0 um³. The image analysis pipeline used to segment the resulting data is built as follows: after stitching all the images together, they run a 3D segmentation algorithm composed of a median noise removal filter, a Frangi vesselness filter detector, adaptive thresholding, and a cluster discard criterion based on signal-to-noise ratio. The resulting binary mask is then skeletonized and turned into a graph [START_REF] Ji | Brain microvasculature has a common topology with local differences in geometry that match metabolic load[END_REF]. While the resolution and precision of the imaging pipeline seem remarkable, it is a bit of an overkill since the minimal capillary radius is around 2 um. The full imaging and processing time is too long, and the imaging is based on gel perfusion which can be prone to inaccuracies due to perfusion quality.

II.5.2 Light sheet microscopy and CLARITY clearing -Di Giovanna et al.,

2018)

For this pipeline, they used the clearing protocol of Clarity, combined with blood vessel lumen staining [START_REF] Di Giovanna | Whole-Brain Vasculature Reconstruction at the Single Capillary Level[END_REF]. Similarly to the previous pipeline, they perfused the mouse with a fluorescent gel comprised of protein to resist the aggressive Clarity lipid washout of the Clarity protocol. Clarity forms a hydrogel-tissue hybrid that retains the proteins and allows selective removal of lipids. The delipidated sample is then incubated in a refractive index matching medium to make it transparent and allow light sheet microscopy acquisition. They imaged the brain at a voxel resolution of 0.65 * 0.65 * 2 um³. Segmentation was performed both automatically based on Markov random field algorithms and manually with commercial software. Not only is Clarity an expensive pipeline, but MRF is usually not enough to propose good segmentation when used alone in a noisy grayscale image. In [START_REF] Di Giovanna | Whole-Brain Vasculature Reconstruction at the Single Capillary Level[END_REF], between 20% and 30% of the voxels segmented by the pipeline are false positive. Moreover, the pipeline does not seem to provide a great variety of analysis tools. Furthermore, Clarity is known to introduce tissue distortion due to the aggressive lipid removal protocol and to be slow.

II.5.3 VesSAP : organic clearing and gel infusions

This pipeline is maybe the closest to ours. It is applied on 3DISCO processed sample images with light sheet microscopy. They segmented the data using a 5-layer Convolutional Neural Network. To circumvent the large amount of data to train the model, they pre-trained the network on a synthetically generated vessel-like dataset and used transfer learning approach so they could refine the training on a small number of manually annotated real vascular data (Todorov et al., 2020). They claimed a 93% accuracy. From our experience, sample quality and signal intensity can vary significantly from one batch to another, making a fully CNN-based pipeline unreliable from one experiment to another without prior refined training on similar data.

II.5.4 Summary of the main pipelines

There are many other pipelines out there using different variations of transparization and analysis protocol.

The detail of those strategies can be found in the table below:

III Network Theory for transportation network

This chapter does not aim to give an exhaustive review of all the graph theory methods but to focus on the general ideas and the ones relevant for studying vasculature. Whether they are used to describe physical structures, social interactions, or financial streams, networks are everywhere in nature and civilization. Spider webs, leaf veins, highways, genomic chemical pathways, and neuronal connections are just a few examples of naturally arising topological wonders. Therefore, many algorithms have been developed to understand these structures, including the vasculature in both physiological and natural transportation networks.

In the case of vasculature, after acquiring and binarizing the vasculature, it is embedded in the form of a graph, which is a mathematical structure of nodes and edges that link nodes together. A conventional formalism represents vascular branch points or vessel bifurcations as nodes, and vessels themselves as edges. The vasculature can be seen as a transportation network. Recent observations tend to describe it as a vastly heterogeneous system with a wide density scale. Since it transports blood, which carries nutrients and oxygen, and considering that these resources are quickly depleted during exchanges between the blood and neural cells, the blood doesn't flow in the capillary network for an extended distance. Instead, it goes in and out of the parenchyma in short loops, and a single red blood cell does not travel across several brain regions in one systolic cycle. A standard intuition in such a physical network is that it can be subdivided or partitioned into homogeneous subpopulations, ideally with low intra-class heterogeneities and high inter-class heterogeneities.

III.1 Overview of common graph analysis strategies

One of the main difficulties in exploring vascular network, is that very little is known about its organization.

However, it is believed that it should be more or less organized into perfusion territories, vein or artery centered, with enough redundancy to keep perfusing any area even after some vessel removal. Because of its numerous applications in urban hierarchy, road tracing, flow monitoring, identifying territories and analyzing the relationships and between them has long been a great deal. Many models and methods emerged to achieve this. Key aspect to identify is :

• redundancy

• territories

• loopiness

Finally, a last way to evaluate it is to compute the eigenvectors of the adjacency matrix. The greatest eigenvalue results in the centrality measure.

One of the main limitations of these concepts is that one centrality will give very different results from another and the choice of one before another differs depending on the application and it only gives relative importance of one node compare to the other with regards of a specific network mechanism assumption (Fig 13).

III.1.2 Networks community

The idea of graph clustering or community partitioning is to divide a large-scale network into subnetworks by maximizing a graph feature that represents the ratio of intra-versus inter-community edges. Here, our aim is to understand the vascular organization. The previously described centrality measures can help as a graph feature in predicting how fluid might be directed through this blood transportation network and identifying the key structures that allow precise distribution and regulation. For instance, vascular networks are subdivided into venous vs arterial domains, or subnetworks irrigated mostly from one of the penetrating arteries, where veins act as "waste" collectors at the frontier [START_REF] Kennel | From whole-organ imaging to in-silico blood flow modeling: A new multi-scale network analysis for revisiting tissue functional anatomy[END_REF] .

Do vascular domains match neuronal substructures? In 2013, Blinder et al. looked at the correlation between fluid flow domains, community partitioning, changes in vascular density, and cortical columns in the barrel cortex (Blinder et al., 2013) . They showed that none of these colocalized with barrel columns, and neither did penetrating arteries and veins, thus strengthening the idea that cerebral vasculature is devoid of partially isolated subnetworks and can be modeled as a regular lattice that homogeneously perfuses the parenchyma [START_REF] Blinder | Topological basis for the robust distribution of blood to rodent neocortex[END_REF][START_REF] Shih | The smallest stroke: Occlusion of one penetrating vessel leads to infarction and a cognitive deficit[END_REF] .

There are many ways to infer communities in graphs, but one of the first problems to address is the definition of a community specific to the type of network, as there is no commonly accepted definition.

From a local standpoint, communities can be built by focusing on internal cohesion criteria that the community finder algorithm will try to maximize. For example, looking at cliques (complete maximal subgraphs with at least three vertices), loop agglomerates, or edge density. From the point of view of vasculature, detecting veins and arterial domains as in [START_REF] Kennel | From whole-organ imaging to in-silico blood flow modeling: A new multi-scale network analysis for revisiting tissue functional anatomy[END_REF] is a relevant process.

From a more global standpoint, the main idea is to compare structured graphs with communities to random graphs, often referred to as Erdos-Renyi graphs, with a certain number of constraints to enable the comparison. In 2004, Newman et al. proposed an implementation of this idea by comparing the graph of interest to a randomized version of it, where the edges were randomly rewired with the constraint of keeping the degree distribution of the original graph (Newman, 2004) . He also introduced the concept of modularity as a measure of the quality of a graph partition. It can be defined as the sum of the probabilities of two adjacent vertices belonging to the same community, summed over all the vertices of the graph [START_REF] Fortunato | Community detection in graphs[END_REF] .

III.1.3 Approaches inspired by vascular biology

Theoretical strategies to analyze graphs often do not match the properties of vascular graphs as defined by their biological mechanisms. This chapter will focus on approaches that may be more faithful to the biological reality of vascular graphs, such as a community finder algorithm called stochastic block modeling. We will compare this to natural communities of vessels, as defined by their arterial or venous domains, as well as some physical models of blood flow through the network.

III.1.3.1 Stochastic Block Model

It is indeed very challenging to describe a vascular structure at large scale, as it contains millions of nodes. In the scope of understanding how vasculature arises and how it rewires with natural plasticity, the interest lies not only in finding its different communities but also in modeling generative mechanisms responsible for network modification.

The most standard Stochastic Block Model (SBM) was initially developed to describe the structures of friendship networks (Holland et al., 1983;[START_REF] Wasserman | Stochastic a posteriori blockmodels: Construction and assessment[END_REF] . The concept of SBM has been adapted to handle more complex networks with a wider degree distribution [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF] , and later

Peixoto covered most of its variants using microcanonical representations to include weighted, oriented graphs, degree correction, and many more (Peixoto, 2014[START_REF] Peixoto | Nonparametric weighted stochastic block models[END_REF][START_REF] Peixoto | Bayesian stochastic blockmodeling[END_REF] .

The basic idea behind SBM is to group nodes with similar connectivity patterns. The original algorithm resulted in a non-convex cluster because high-degree nodes tended to cluster together. Therefore, a degree correction term representing the degree of each node was introduced in the likelihood expression to allow for heterogeneities within each cluster [START_REF] Funke | Stochastic block models: A comparison of variants and inference methods[END_REF] . Roughly, the implementation is as follows: each node is randomly assigned a cluster and then randomly moved to another cluster. The change is accepted or rejected according to a probability law. The clusters are also randomly merged, and the merge is accepted or rejected based on some probability. To e v a l u a t e t h e q u a l i t y o f a p a r t i t i o n , t h e c o n c e p t o f m o d u l a r i t y w a s i n t r o d u c e d . I t m e a s u r e s t h e e x t e n t t o w h i c h edge density within a cluster is higher than between clusters. However, it has been shown that in lattice and regular graphs, even if they are random graphs, the modularity remains quite high due to their intrinsic lattice topology. For instance, in a cubic graph like the vasculature, modularity is likely to fall within the interval [0.66- [START_REF] Mcdiarmid | Modularity in random regular graphs and lattices[END_REF] . Modularity is a useful tool for inferring vascular groups and characterizing the brain vascular network. However, it's important to remember that any spatially convex partition of this network will be highly modular, even if it doesn't correspond to real substructure. Maximizing the modularity of the vascular graph becomes essential to ensure that the resulting partition is meaningful.

0.88] (McDiarmid

III.1.3.2 Mecafluidic simulation/quantification

Since the vasculature is a transportation network flow velocity and pressure gradient for the extracted fraction are important parameters to evaluate. Flow modeling has become increasingly important as a research field owing to the availability of high-quality 3D segmentations of large vascular networks. They can be used as predictive models to investigate different scenarios of pathological evolutions for instance. If the blood is considered a pure plasma, with a given geometry and boundary conditions the flow, pressure, stress, and

This table summarizes the main variations of the Stochastic Block Model algorithm. Those variants may consider the degree distribution of vertices in each cluster to avoid clustering high degree vertices which is a shortcoming of the standard SBM method. In other variants, the user can consider edge weigh to control the clustering. Another interesting feature is the possibility to set the number of clusters. This is particularly important to benchmark for modularity, or if many modules are needed. Hierarchical SBM perform standard SBM at different scales recursively. All these techniques have their pros and cons. In the case of this work, the standard SBM was used to avoid biasing the module search.

other parameters are the solution to the Navier Stokes equation which is a set of partial differential equations representing the balance between, pressure, viscous and gravitational forces. To solve the system, a common method is to use simplification.

In the equation below (Eq 1):

u represents the fluid velocity, ρ is the volumic mass, p is the fluid pressure, μ is the fluid viscosity, and f the gravity.

Under certain conditions some of the terms of the Navier Stokes equation can be neglected. One of the parameters governing this is the Reynolds number. It represents the ratio between inertial and viscous forces.

High values of Reynold number indicate a turbulent flow, and lower values indicate laminar flow. For liquid such as blood, whose viscosity is only 2 to 3 time superior to water viscosity it is the fact that it is circulating through extremely narrow tubes (a capillary vessels diameter is in average 2 to 5 um) that enable us to consider it has a low Reynolds number (Eq 2), thus allowing to drop some of the Navier Stokes equation terms and thus simplified the resolution, especially if it combines with the steady state assumption [START_REF] Beier | Vascular Hemodynamics with Computational Modeling and Experimental Studies[END_REF] .

Boundary conditions are also essential in the correctness of the model. The most common is to choose the inlet or outlet flow pressure and velocity profiles (Eq 3). This is accessible thanks to in vivo measurement done in surface arteries or veins that are accessible [START_REF] Beier | Vascular Hemodynamics with Computational Modeling and Experimental Studies[END_REF] .

This allows one to compute the pressure and flow velocity at any point of the vascular network and to predict the effect of removal or clogging of any vessels on the general blood circulation [START_REF] Pries | Blood Flow in Microvascular Networks[END_REF]) (Hirsch , 2012) , or the effect of the lattice development on the flow velocity. For instance, in (Goddet, n.d.) , it is showed that in an electrical mesh, the tighter the mesh, the slower its intensity. With the previously described simplification, the Navier Stokes equation can be seen as a set of linear equations, similar to the Ohm law where flow is equivalent to intensity, and pressure to potential. Immediately, it means that a densely branched vascular network will see a slow flow. This was observed for the vasculature in [START_REF] Mirzapour-Shafiyi | Numerical evaluation reveals the effect of branching morphology on vessel transport properties during angiogenesis[END_REF] , and in [START_REF] Hirsch | Topology and hemodynamics of the cortical cerebrovascular system[END_REF] .

Formulation of the flow computed by Poiseuille's law, and equivalent resistance empirically derivated by Pries and Secomb in 2005 (Eq 4):

In [START_REF] Secomb | Blood Flow in the Microcirculation[END_REF] the author also discusses the network approach to studying blood flow in the microcirculation. The network approach involves considering the microcirculatory system as a network of interconnected vessels, ranging from larger arterioles to smaller capillaries and venules. This approach focuses on analyzing the structure and organization of the microvascular network, including parameters such as vessel diameters, branching patterns, and connectivity. By adopting a network perspective, researchers can investigate how the geometry and topology of the microvascular network influence blood flow distribution and overall hemodynamics. This approach allows for the examination of flow heterogeneity, pressure gradients, and flow redistribution within the network. Additionally, the network approach can be used to study the effects of vascular remodeling, such as the formation of collateral vessels, on blood flow dynamics. Through mathematical modeling and simulation techniques, the network approach provides a quantitative framework to understand the complex behavior of blood flow in the microcirculation. It offers insights into phenomena such as flow regulation, oxygen delivery, and the impact of diseases or interventions on microvascular function.

By combining experimental data and theoretical analyses, the network approach contributes to a deeper understanding of the microvascular system's role in physiological processes and pathological conditions. In the same paper, the author discusses the hematocrit approach to study blood flow dynamics in the microcirculation. The hematocrit approach is based on the understanding that blood consists of plasma and red blood cells (RBCs), with RBCs occupying a certain volume known as hematocrit. This approach considers the effects of RBCs on blood flow by accounting for their rheological properties, such as their concentration and deformability. By incorporating these factors, the hematocrit approach enables the evaluation of the distribution of RBCs within blood vessels and their impact on flow resistance and velocity profiles. Additionally, this approach allows for the assessment of oxygen transport and exchange, as well as the influence of hematocrit on various physiological processes. Overall, the hematocrit approach provides valuable insights into the complex interplay between blood flow dynamics and the cellular components of blood in the microcirculation.

Another application is to understand the organization of the cerebral blood flow within the different brain regions. In [START_REF] Schmid | Depth-dependent flow and pressure characteristics in cortical microvascular networks[END_REF] , they used a simplified version of the Navier Stokes equation to model and compute the flow within each vessel, and used several models of boundary condition, by embedding their chunk of vasculature inside a larger artificial lattice network. Within the cortex flow seems to show laminar behavior through each layers, mostly driven by pressure drops along the way.

III.1.3.3 Single vessels statistics: tortuosity, orientation

Initially, the vascular network could be thought of as isotropic and homogeneous, which would greatly facilitate the design of a realistic model of the capillary bed [START_REF] Peyrounette | Multiscale modelling of blood flow in cerebral microcirculation: Details at capillary scale control accuracy at the level of the cortex[END_REF]. In [START_REF] Ji | Brain microvasculature has a common topology with local differences in geometry that match metabolic load[END_REF]) , the regional anisotropy of the network was quantified using Monte Carlo sampling. It was found that the heterogeneity was particularly high in the Superior colliculus, pons, and corpus callosum. These types of quantifications, coupled with flow information, for instance, would help understand the relationship between flow control and directionality, vascular and neuronal density, neuronal demand, and capillary orientation (Fig 15).

Figure 15: from Pries et al 1990 -Schematic drawings showing possible reasons for high resistance in small microvessels. In flow calculation the geometry of the vessel itself (namely the radius) is not the only parameter to consider. The tortuosity or irregularities of the vessels can lead to slower blood flow and altered nutrient diffusion. 90-degree intersection, a locally constricted vessel cross section, can lead to an accumulation of cells slowing down the blood flow. Even for models that don't take these parameters into account, it is important to keep in mind that the dynamic of the blood circulation and its directionality can be strongly impacted by the small irregularities, and therefore vessels have to be individually characterized instead of being considered as an homogeneous mass.

III.1.3.4 Diffusion and blood supply

Since the vascular network fills the whole parenchyma to efficiently distribute oxygen and nutrients, a relevant measure is the distribution of minimal distances between any point of the brain tissue and the closest vessel lumen. [START_REF] Ji | Brain microvasculature has a common topology with local differences in geometry that match metabolic load[END_REF] found that the Inferior Colliculus, a region known to have a high vessel length density, display the smallest distances between tissue and vessels in average. In the opposite site of the spectrum is the corpus callosum, a myelinated tract.

They also proposed to compare the relationships between vessel density and tissue to vessel distance in different lattice models, and showed that by modifying its topology, the capillary network can reduce the tissue to vessel distance while keeping its vessel density unchanged. This set the basis to explore the effect of single capillary clogging or removal on the tissue access to oxygen and nutrients, and thus to maybe investigate the link between neurodegenerative diseases and vascular reduction.

The nutrients flow through the vessels to the brain via diffusion through the vessel membrane. The idea is the see the relation between vessel density, and tissue to vessel distance, and the diffusion of these in the tissue [START_REF] Weber | The microvascular system of the striate and extrastriate visual cortex of the macaque[END_REF] . This can be described using the concept of extracted fraction. If the solute is permeable, it crosses the BBB at a certain rate and diffuses in the tissue (Fig 16). The diffusion equation is driven by the flow, which is itself driven by the neural activity, the geometry of vessels, and the diffusion coefficient. Using the flow equations previously described, one can compute its value on the entire brain, and thus could infer useful information about the correlation with neural activity (Q. R. [START_REF] Smith | A review of blood-brain barrier transport techniques[END_REF] .

Figure 16: nutrients and oxygens are distributed to the brain via diffusion through the membrane of the blood vessels. The diffusion rate can be modeled as follows: nutrient concentration [C] follows a law of exponential decay. The uptake is (Ca-Cv) * F , where Cv is the concentration downstream and Ca the upstream concentration. F is the flow and can be computed as the time it takes for a volume V to run through the vessel. The uptake also depends on the membrane surface the blood is in contact with. To summarize, uptake increases in large vessel with low flow. Since arteries and veins have higher flow rates, most of the diffusion occurs in capillaries where the flow is much slower. Computational models show that the more branches the capillary bed has, the slowest is the flow, and therefore the highest the diffusion. Therefore, vessel density is a good indicator of areas with large nutrient intake.

III.2 Generation of artificial vascular networks in silico

Applying algorithms to vascular graphs is not straightforward. In fact, biological vascular graphs reach very large size, namely up to 3 to 4 million nodes and edges and require highly effective algorithms that can scale.

In addition, vascular graphs present the following characteristics: they are sparse, 3-regular graphs, making betweenness centrality, and hub searching quite irrelevant.

III.2.1 Degree 3 rules

Because of the way the cerebral vasculature matures (iterative rounds of sprouting), it only contains degree 3 nodes [START_REF] Kur | Temporal modulation of collective cell behavior controls vascular network topology[END_REF] . According to [START_REF] Fleury | Diffusion limited aggregation from shear stress as a simple model of vasculogenesis[END_REF], 2000;[START_REF] Nguyen | Dynamics of vascular branching morphogenesis: The effect of blood and tissue flow[END_REF] , vascular networks develop through the iterative invasion of penetrating veins and arteries into the already formed capillary bed.

This leads to the model of vasculature as fractal trees linked together via a homogeneous 3-regular capillary mesh. The first convincing 2D model for vasculature was based on Voronoi tessellation, which served as a scaffold for the capillary mesh, thus respecting the degree 3 constraint of vascular networks [START_REF] Lorthois | Fractal analysis of vascular networks: Insights from morphogenesis[END_REF] .

These models aim to describe the vascular network from a physical point of view and reproduce its physiological properties. However, they do not attempt to replicate how the networks are set up in real life. For a growth-realistic model, Vanessa Coelho et al [START_REF] Coelho-Santos | Postnatal development of cerebrovascular structure and the neurogliovascular unit[END_REF] showed that inside the basic capillary plexus, some vessels mature early as arteries and veins. Later, the capillary network densifies through the sprouting of new vessels branching from venous domains to reconnect to nearby arterial vessels. This process is controlled via a VEGF/DLL4 gene network [START_REF] Kur | Temporal modulation of collective cell behavior controls vascular network topology[END_REF] . In addition, [START_REF] Mada | Pattern formation of vascular network in a mathematical model of angiogenesis[END_REF] proposed a fractal model to mimic the retinal vasculature pattern using a simple set of differential equations that describe the dynamic growth of the network, taking into account the cell division rate, distribution of branching angles, and vessel connectivity.

III.2.2 Topological specificity of the vascular networks

Until recently, not many studies focused on cerebral vascular topology, and the capillary mesh was often simplified as a regular mesh or only veins and arteries were used [START_REF] Lorthois | Fractal analysis of vascular networks: Insights from morphogenesis[END_REF][START_REF] Peyrounette | Multiscale modelling of blood flow in cerebral microcirculation: Details at capillary scale control accuracy at the level of the cortex[END_REF].

in addition to being 3-regular, the network of capillary vessels has the property of being space-filling [START_REF] Daniels | The Nucleotide Sensor ZBP1 and Kinase RIPK3 Induce the Enzyme IRG1 to Promote an Song et al[END_REF], which is a characteristic shared with 3D Voronoi diagrams. However, 3D Voronoi graphs show highdegree nodes (5 and more neighboring nodes). In [START_REF] Goirand | Network community structure and resilience to localized damage: Application to brain microcirculation[END_REF] , they used random merge, pruning, and addition of vessels to fit Voronoi graphs into the vascular constraints and reproduce the physical, topological, and transportation properties of the brain capillary network. They also measured the modularity of such graphs and compared it with that of the mouse brain cortical areas. Indeed, such partitioning property has already been shown in vascular networks (Blinder et al., 2013) , with a modularity value of 0.946. With

Voronoi, they were able to reach 0.827, which is much closer to the average modularity of a random 3-regular graph (~0.661) [START_REF] Goirand | Network community structure and resilience to localized damage: Application to brain microcirculation[END_REF] . This result confirms that the vasculature forms spatially segregated clusters of vessels due to its space-filling property. The cerebral vasculature is a very dense network of capillaries, veins, and arteries, strongly intertwined with other cerebral functional networks such as neuronal systems, lymphatic systems, glial systems, and so on. It is known to be an active component in numerous brain diseases, but reconstructing vasculature at the wholebrain scale and quantifying interindividual variations from intergroup variations is extremely challenging.

IV.1. Summary of

Combining tissue clearing, immunolabeling, and novel computational techniques, we propose a pipeline to image, process, and analyze terabytes of light sheet microscopy multichannel data.

The pipeline is divided into 3 main parts: on the biological side, podocalyxin and CD31 were used to stain the whole vasculature, and Acta2 was used to stain the arteries and veins. A novel Wobbly stitcher was developed to fix the shifts introduced by the light sheet engine axis, and multiple filters were combined in order to capture every capillaries and veins despite the wide range of different signal intensities and scales.

Last but not least, we use deep neural networks to correct for artifacts introduced by the immunolabeling nature and prevent topological aberrations.

We built a dataset of over 20 animals' brains of millions of vertices vascular graphs, with labeled veins and arteries. Using this dataset, we characterized the heterogeneities and organizational patterns of the brain vasculature across different regions, proposed vascular topology-based classifications, and analyzed the rearrangement of vessels in pathological diseases such as strokes.

We combined graph theory algorithms and physical modeling to characterize the modularity, the arteriovenous vascular domains, and blood circulation paths in control and pathological conditions. We used a congenitally deaf mice model to quantify the effect of immediate and long-range effects of neuronal activity loss on the vascular structure.

In Brief

TubeMap, a pipeline for characterizing brain-wide variations in vascular topology, including vessel network structure, geometry, and arterio-venous labeling, was developed and used to investigate vascular remodeling in models of ischemic stroke and sensory loss.

INTRODUCTION

Blood flows within the brain through an intricate vascular network of arteries, capillaries, and veins. The high metabolic needs of the brain and its low capacity to store energy require that neuronal and vascular functions are intimately linked [START_REF] Girouard | Neurovascular coupling in the normal brain and in hypertension, stroke, and Alzheimer disease[END_REF]. Minimal perturbations of the blood flow impact neuronal function and can compromise neuron survival [START_REF] Zlokovic | The blood-brain barrier in health and chronic neurodegenerative disorders[END_REF]. Conversely, neuronal activity modulates hemodynamics [START_REF] Ances | Coupling of changes in cerebral blood flow with neural activity: what must initially dip must come back up[END_REF]. Thus, many neurological disorders have a vascular component, and most vascular alterations have consequences on brain functions [START_REF] Iadecola | The pathobiology of vascular dementia[END_REF][START_REF] Kisler | Cerebral blood flow regulation and neurovascular dysfunction in Alzheimer disease[END_REF][START_REF] Zlokovic | Neurovascular pathways to neurodegeneration in Alzheimer's disease and other disorders[END_REF].

The cerebral vasculature is challenging to study as it spans several spatial scales (from micron-sized capillaries tiling the brain to vessels extending over several millimeters) and forms a highly complex network [START_REF] Begley | Structural and functional aspects of the blood-brain barrier[END_REF]i n which exchanges with the brain are tightly controlled by the blood-brain barrier [START_REF] Zlokovic | The blood-brain barrier in health and chronic neurodegenerative disorders[END_REF]. 2-photon microscopy combined with cranial windows has been an essential tool to study the dynamics of blood flow and metabolism in vivo [START_REF] Ma ¨chler | In Vivo Evidence for a Lactate Gradient from Astrocytes to Neurons[END_REF][START_REF] Roche | In vivo imaging with a water immersion objective affects brain temperature, blood flow and oxygenation[END_REF]. However, the impact of the vascular topology on the function of neural circuits is often overlooked because of the difficulty to extract information on the brain-wide vascular organization. Addressing these questions requires systematic comparisons of the fine structure of the vasculature in large volumes.

A comprehensive description of the organization and remodeling of the vasculature also calls for methods allowing the identification of the arterial, capillary, and venous components. Large-scale cerebral vascular reconstructions have been carried out to this day predominantly with dye filling strategies, in which a vessel-filling solution is perfused intracardially (Blinder et al., 2013;[START_REF] Kleinfeld | Large-Scale Automated Histology in the Pursuit of Connectomes[END_REF]Quintana et al., 2019;[START_REF] Tsai | Correlations of neuronal and microvascular densities in murine cortex revealed by direct counting and colocalization of nuclei and vessels[END_REF]. More recently, such injections have been combined with light sheet microscopy in optically cleared samples for large-scale reconstructions [START_REF] Breckwoldt | Correlated magnetic resonance imaging and ultramicroscopy (MR-UM) is a tool kit to assess the dynamics of glioma angiogenesis[END_REF][START_REF] Di Giovanna | Whole-Brain Vasculature Reconstruction at the Single Capillary Level[END_REF][START_REF] Lagerweij | Optical clearing and fluorescence deep-tissue imaging for 3D quantitative analysis of the brain tumor microenvironment[END_REF][START_REF] Liebmann | Three-Dimensional Study of Alzheimer's Disease Hallmarks Using the iDISCO Clearing Method[END_REF][START_REF] Lugo-Hernandez | 3D visualization and quantification of microvessels in the whole ischemic mouse brain using solvent-based clearing and light sheet microscopy[END_REF][START_REF] Todorov | Automated analysis of whole brain vasculature using machine learning[END_REF][START_REF] Zhang | CLARITY for High-resolution Imaging and Quantification of Vasculature in the Whole Mouse Brain[END_REF]. While vessel-filling allows for a bright labeling that facilitates imaging and downstream analysis, the nature of the reconstructed vessels is unknown, and can, to this day, only be retrieved manually on highly resolved datasets (Blinder et al., 2013). Additionally, reproducibility is challenging with intracardiac perfusions (Quintana et al., 2019). Alternative labeling methods have been used, such as tail-vein lectin injections (d'Esposito et al., 2018), or the use of mouse genetic lines with labeled vessel walls [START_REF] Jing | Tissue clearing of both hard and soft tissue organs with the PEGASOS method[END_REF]. However, because these labels are faint, discontinuous, and leave the vessel centers Cell 180,12345678910111213141516February 20,2020 ª 2020 Elsevier Inc (legend continued on next page) 2 Cell 180,12345678910111213141516February 20,2020 unlabeled, they pose an unsolved challenge for image processing at high resolution and large scales.

We explored the possibility of using immunolabeling and tissue-clearing to construct graphs of vascular networks. This strategy, in principle, offers several advantages. The structural labeling of vessel walls does not require intracardiac perfusion and would therefore allow the labeling of organs or embryos that are difficult or impossible to perfuse. Moreover, immunolabeling would enable the combination of a generic vessel signal with markers specific for each vessel type. So far, it has been difficult to balance labeling specificity versus epitope density using immunolabeling against the vessel walls in intact adult organs. In addition, there is also a lack of tools designed to process data from 3D fluorescent microscopy on terabytesized image stacks. Light sheet imaging and tiling artifacts lead to frequent duplication or fusions of vessels that greatly alter network topologies and compromise large-scale statistical comparisons of vascular datasets.

Here, we report a method to generate labeled (arteries, capillaries, and veins) and region-annotated vascular graphs of the mouse brain. It relies on an immunolabeling cocktail that produces a bright labeling of the vascular endothelial and smooth muscle cells in unsectioned brains using the iDISCO+ technique (Renier et al., 2016). To construct vascular graphs, we designed a complete processing pipeline, TubeMap, which addresses issues of vessel duplication and other imaging artifacts as well as the multi-scale nature and hollow tube structure of the vasculature network. TubeMap is optimized to work on large terabyte datasets and integrates functions to handle, visualize, and analyze the reconstructed vascular networks made of several millions of vessels.

Using this novel pipeline, we measured the regional variations of the mouse cerebral vascular topology. We found that while the density of the capillaries can vary up to an order of magnitude between regions, their distance to the nearest arterial or venous endpoints are similar throughout the brain. We used clustering methods to compare the organization of neocortical regions. Finally, we analyzed the remodeling of the vascular network in two different pathological conditions: congenital sensory deprivation of the auditory system and brain ischemia in adult animals. These models showed two distinct modes of vascular remodeling and plasticity, with cross-modal long-range changes of the vascular density in congenitally deaf mice and region-independent reorientation of capillaries in stroked brains.

RESULTS

Labeling and Imaging of the Vasculature Complete labeling of the brain endothelial cells with antibodies is challenging, due to the density of epitopes. We verified that the iDISCO+ processing time and pre-treatments did not affect the integrity of the vasculature (Figure S1A). To enhance the signal of the vessels immunolabeling, we sought to design a ''super-polyclonal'' immunoglobulin cocktail and tested a panel of antibodies targeting proteins broadly expressed by endothelial cells: CD31 [START_REF] Newman | PECAM-1 (CD31) cloning and relation to adhesion molecules of the immunoglobulin gene superfamily[END_REF], collagen IV, and podocalyxin [START_REF] Horvat | Endothelial cell membranes contain podocalyxin-the major sialoprotein of visceral glomerular epithelial cells[END_REF][START_REF] Testa | Ubiquitous yet distinct expression of podocalyxin on vascular surfaces in normal and tumor tissues in the rat[END_REF]( Figure 1A). None of the anti-CD31 antibodies tested yielded sufficient signal. Anti-podocalyxin labeling improved the signal deep in the brain, but still with a low signal-to-noise ratio at the center when used on its own (Figures 1C andS1B). However, combining both markers labeled the vasculature with high signal-to-noise ratios throughout the brain (Figure 1C). To increase the flexibility of the method, we tested different host combinations of these markers (Figures S1B-S1D; Methods S1). We verified that the complete vasculature was labeled with this strategy, by staining for CD31 on sections after an anti-podocalyxin whole-mount labeling (Figure S1E).

We complemented this labeling with additional markers specific to arteries. We explored the use of smooth muscle actin (Acta2) and transgelin (Sm22), both known to have high expression levels within arterial walls and weaker in veins [START_REF] Brunet | Netrin-1 controls sympathetic arterial innervation[END_REF][START_REF] Hill | Regional Blood Flow in the Normal and Ischemic Brain Is Controlled by Arteriolar Smooth Muscle Cell Contractility and Not by Capillary Pericytes[END_REF][START_REF] Vanlandewijck | A molecular atlas of cell types and zonation in the brain vasculature[END_REF]. iDISCO+ labeling of Acta2 and Sm22 yielded, as expected, a dense arterial signal, while surface and large penetrating veins had a weak and patchy distribution (Figure S2A; Video S1). Validation of the volume Acta2 labeling on sections confirmed that all arteries were labeled but smaller, higher order branches were occasionally only visible from the section staining (Figure S2B). We finally tested the possibility to discriminate veins by screening for venous-enriched markers [START_REF] Vanlandewijck | A molecular atlas of cell types and zonation in the brain vasculature[END_REF] and found that von Willebrand factor (vWF) immunolabeling was compatible with volume staining and enabled veins/arteries discrimination on raw data (Figure S2C).

Vascular imaging requires high resolutions, and we strived to find a compromise between resolution, data handling, and acquisition speed using a commercial light sheet microscope. To maximize the z resolution, we optimized the sheet curvature and tiling layout (Figure 1B). The apparent diameters of Cell 180, 1-16, February 20, 2020 3 capillaries on raw data were between 2.6 and 2.7 px in x or y (about 4.2 mm) and between 4.0 px in z at the center of the tile to 5.8 px at the level of the blended tiles overlap, the less resolved portions of the data (Figures S3A-S3C).

Stitching the mosaic accurately is necessary to prevent capillary duplications. Commonly used stitchers, such as Terastitcher [START_REF] Bria | TeraStitcher -a tool for fast automatic 3Dstitching of teravoxel-sized microscopy images[END_REF] or Imaris Stitcher produced many tiling errors on our data (about 10% of the borders where misplaced) (Figures S4A andS4B). Although BigStitcher (Ho ¨rl et al., 2019) improved on tile placement, we noticed that mechanical jitters of stage movements still prevented accurate alignments throughout the stacks (Figure 1F). To solve these problems, we developed WobblyStitcher, a tool for robust and non-rigid stitching of terabyte datasets (Figures 1G and1H) (see STAR Methods for details on the stitcher and computing resources). WobblyStitcher produced accurate non-rigid placement of the data and was robust in placing challenging planes at the edge of the sample, virtually eliminating object duplications in the dataset and producing continuous images of vessels (Figures 1I, S4C, andS4D).

Vasculature Graph Construction

Despite a high signal-to-noise ratio (Figure S3A), we found that a processing pipeline based on a single filter could not accurately segment vessels in our data due to glow and shadows artifacts, but also because of the widely varying sizes and signal intensities of vessels. To tackle the multi-scale nature of our data, we implemented a collection of filters arranged in parallel paths (''multi-paths'') designed to process complex fluorescent microscopy images of very large sizes (TB) in a few hours on a standalone workstation (Figures 2A and2B). Lightsheet imaging artifacts of streaking shadows that resemble vessels are first removed with 2 crossed 3D rank filters (Figure 2C). In a second step, the multi-path pipeline is used to binarize the vessels. Each path captures and handles different aspects of the vessel signal (Figures 2D-2G).

This first part of the pipeline can binarize different types of vascular datasets generated from gel filling or wall labeling methods. However, the labeling of endothelial and arterial walls produces images of both solid (capillaries) and hollow tubes (veins and arteries). The segmentation of hollow vessels has not been attempted yet on very large data [START_REF] Bates | Segmentation of Vasculature From Fluorescently Labeled Endothelial Cells in Multi-Photon Microscopy Images[END_REF]. We used a machine-learning approach and designed a deep convolutional neural network (CNN) based on an architecture optimized for vessel segmentation on large 3D datasets to fill empty vessels [START_REF] Livne | A U-Net Deep Learning Framework for High Performance Vessel Segmentation in Patients With Cerebrovascular Disease[END_REF]. Our CNN was trained using whole-brain datasets containing only filled vessels and imaged using the binarization pipeline (Figures 2H, S5A, and S5B; STAR Methods). Tested on an artificial dataset made from a real graph of vessels with radii ranging from 2 to 25 pixels (Figures S5C andS5D), the network had a false positive rate of 0.1% ± 0.06% pixels. We applied this CNN to both the vessels and arterial binary masks to obtain filled tubes throughout the datasets (Figures 2H, S5E, andS5F).

The binary image is then used to construct the graph of the vasculature network (Figure 2I). We built a highly optimized 3D skeletonization algorithm based on a directional parallel thinning method and discrete topology (Figure 2I). Our skeletoniza-tion step traces the center lines of vessels in a typical dataset of 10 billion foreground pixels in about 1 h.

Graphs were then created from the centerline image and augmented by information about vessel geometry, identity, and brain annotation (Figures 2I andS6; STAR Methods). To determine vessels identity, seed points were placed on clearly identified large veins, based on their radii and low Acta2 expression. For arteries, seed points were placed on vessels with high Acta2 expression. Finally, the labels were expanded by tracing from these seed points down to capillary radii (Figures S6B-S6E).

To validate our pipeline (Figure 3A), we performed manual segmentation from the raw data in random cubes extracted from 3 locations in a whole brain scan, representative of the data heterogeneity (Figure S5G). After the manual segmentation, the whole scan was processed with TubeMap and the graph constructed. We found that 4.0% ± 1.7% of manually annotated vessel branches were missed by TubeMap, and 2.9% ± 1.8% of the branches present in the TubeMap graph had not been manually annotated.

Analysis of Vascular Graphs

We first constructed vascular graphs of brains from three perfused C57Bl6 3-month-old mice labeled for podocalyxin, CD31, and Acta2. We plotted 3D renders of these graphs to verify the consistency of our unsupervised reconstructions from identifiable hallmarks of the rodent's cerebrovasculature (Figures 3B-3F). In the cortex, alternating penetrating arteries and veins as well as the laminar heterogeneity of capillary density were clearly visible, with higher densities of vessels in the isocortex than in the allocortex (Figure 3B) [START_REF] Michaloudi | Areal and laminar variations in the vascularity of the visual, auditory, and entorhinal cortices of the developing rat brain[END_REF]Schmid et al., 2019). Known laminar or radial patterns of capillary densities were visible in the hippocampus (Figure 3B), brainstem, midbrain (Figure 3D), or cerebellum (Figure 3E). We verified the quality of the vessel identity labeling from the expected organization of veins and arteries in the hippocampus, cerebellum, and spinal cord (Figure S7). In the hippocampus, alternating arcs of veins and arteries were observed (Figure 3C). In the cerebellum, arteries were detected primarily at the pial surface of the lobules and veins mostly in the white matter (Figure 3E). In the spinal cord, arteries where visible as arbors fanned out from the anterior and the posterior lateral spinal arteries (Figure S7B; Video S2). Vessels belonging to white matter tracts (Figure 3F, in white) were longer and less dense than their gray matter counterpart.

We then turned to a brain-wide quantitative analysis of the vascular network. The total length of the vasculature added up to 144 ± 2 m per hemisphere (n = 3). Using the Allen Brain Atlas annotation, we further resolved the vessel length density across different brain areas (Figure 4A) with the cortex contributing the largest part. Vessel length density varied widely at a sub-region level: for example, all cortical areas showed repetitive variations in the vessel length density across cortical layers (Figure 4A,inset).

We then included topological aspects into our analysis. We identified all branch points at which vessels split or join. A typical graph consisted of 3.2 million branch points and 4.4 million vessel segments per hemisphere. The brain-wide mean vessel branch point density was 6,400 per cubic millimeter (Figure 4B). However, many brain regions deviate by a wide margin from this mean. The cortex, hippocampus, and striatum had densities below average in most regions with notable exceptions, for instance in the layer 4 of primary sensory regions. The densest nucleus of the brain was the median mammillary (Figures 4A and4B; Table S1). Other notable nuclei known for their high densities are seen among sensory relays or neurosecretory regions and are indicated in Figures 4B and4C.

Notably, the branch point density closely followed the vessel length density (Figures 4A versus 4B). Indeed, the number of vessel segments strongly correlated with the number of branch points in different brain regions with a ratio of 1.78 between them on average (r 2 = 0.94) (Figure 4C). If the vasculature network was a perfect tree in which each vessel would split into two vessels at each branch point (or conversely two vessels would join into one), this ratio would be expected to be precisely 3/2. The measured ratio is larger, indicating a level of redundancy in the vasculature network due the presence of either higher order branch points or loops. We therefore looked at the distribution of 4-points motifs across the brain to search for loops and found that although 85% of the network is made of simple forks, the remaining 15% contain loops indicative of the redundancy of the capillary bed (Figure 4E).

These results confirm that the capillary density of is highly heterogeneous throughout the brain [START_REF] Craigie | The architecture of the cerebral capillary bed[END_REF][START_REF] Michaloudi | Areal and laminar variations in the vascularity of the visual, auditory, and entorhinal cortices of the developing rat brain[END_REF]Schmid et al., 2019). Variations in capillary density are not correlated with neuronal density but could be related to oxidative metabolism at synapses [START_REF] Keller | Vascularization of cytochrome oxidase-rich blobs in the primary visual cortex of squirrel and macaque monkeys[END_REF]. Apart from neurosecretory or homeostat-regulating nuclei, high capillary densities could denote higher levels of oxidative metabolism in synaptic terminals-rich regions that require high concentrations of ATP [START_REF] Ashrafi | Glucose metabolism in nerve terminals[END_REF], as in sensory pathways.

As vessel length and vessel radius determine the resistance for a fluid to pass through a tube, we further measured length and radii of the vessel segments between two branch points. The length distribution in the cortex peaks at 20 mm, in agreement with previous studies (Blinder et al., 2013). However, the cortical distribution had longer vessel segments than other regions (Figure 4F). The joint distribution of length and radii across brain regions was comparable with previous studies in the cortex (Blinder et al., 2013), but again their distributions are slightly different in the other brain areas (Figure 4G).

We then included the arterial and venous annotations into our analysis. We asked whether the spatial variations of branch point and length densities across the brain affect the distances of capillaries to the next artery and vein. For this, we measured the length and the number of branch points along the shortest path from each vessel to the nearest artery or vein (Figure 4H). Interestingly, on average, the distance, as well as the branch order to the next vein, is larger than to the closest artery, suggesting that the distribution of blood requires a finer control than the collection (Figures 4I-4K). On average, a capillary is 5 or 7 branch points away from an artery or a vein, respectively (Figure 4I), and the mean shortest path between an artery and a vein that passes through a given vessel segment is 120 mm long (Figure 4J). We then focused on the distribution of arterial domains in the cortex. We isolated individual cortical arteries (Figure 4L) and estimated their capillary domain based on the previously calculated shortest distances (Figure 4M). We then counted the vessels in each arterial domain. The distribution of domain sizes peaked between 200 to 1,000 vessels per domain, but very large domains made of over 5,000 vessels were also found (see example in Figure 4M).

We finally chose to investigate in more detail the regional variations in the organization of the isocortex by extracting and comparing a selected set of features: we first extracted branch point densities in each layer, which were higher in the sensory cortices, especially in layer 4 (Figure 5A). We also extracted the position of arterial branches and end points (Figure 5B). Next, we measured the orientation of vessel branches, by considering the vector between two branching vertices and determined whether it was predominantly planar or radially oriented. We noticed that in all sensory regions, the upper layers had higher proportions of radially oriented vessels (Figure 5C). Conversely, integrative and motor regions had an equal mix of radial and planar oriented vessels between upper and deeper layers (Figures 5C and5D). To combine the data from all the extracted features, we measured the differences between the distribution of these features across layers averaged from 3 brains. This separated two groups of regions: a block of sensory areas and another of integrative/motor areas (Figures 5D and5E). This finding suggests that the type of computation performed across cortical layers strongly influences the organization of the vascular topology in the cortex.

Mapping Vascular Network Plasticity

The brain vasculature is a plastic network that can incur largescale rearrangements following various types of insults, such as (B-F) Slices through a vascular brain graph, obtained from a C57bl6 3-month-old female mouse, plotted with ClearMap 2.0. Arteries are indicated in red, putative veins in blue, and capillaries are colored according to their regional annotation. (B) 250-mm coronal and sagittal slices of the graph, showing the general organization of the cerebral vasculature and its regional variations. The inset shows a detail at the level of the hippocampus, with an arterial arc visible. (C) 750-mm coronal and sagittal slices taken at the same levels, showing only labeled veins and arteries. Insets show sagittal and coronal views of the hippocampus, showing alternating arterial and venous arcs. (D) 250-mm projection at the level of the brainstem midline, showing the radial organization of veins and arteries in the hindbrain and superior colliculus. (E) Detail of a 1-mm midline sagittal slice at the level of the cerebellum: veins are predominantly located in the white matter and arteries at the pial surface of the lobules. (F) Detail of the anterior commissure at a midline level (sagittal section). The white annotated vessels (fiber tracts) are longer and less dense that the gray matter vessels (red, thalamus; light blue, striatum). Abbreviations: SO, stratum oriens, SP, stratum pyramidale, SLM, stratum lacunosum moleculare, SM, stratum moleculare, SG, stratum granulosum, H,: hylus. Scale bars represent 1 mm, except in (F) (250 mm).

Cell 180,12345678910111213141516February 20,2020 7 Please cite this article in press as: Kirst et al., Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature, Cell (2020), https:// doi.org/10.1016/j.cell.2020.01.028 traumatic brain injury [START_REF] Cao | Three-dimensional imaging of microvasculature in the rat spinal cord following injury[END_REF][START_REF] Sosa | Low-level blast exposure disrupts gliovascular and neurovascular connections and induces a chronic vascular pathology in rat brain[END_REF], stroke [START_REF] Liu | Vascular remodeling after ischemic stroke: mechanisms and therapeutic potentials[END_REF], or epilepsy [START_REF] Arango-Lievano | Topographic Reorganization of Cerebrovascular Mural Cells under Seizure Conditions[END_REF]. It can also be remodeled by changes in neuronal activity levels, albeit mostly during embryonic and early post-natal development (Harb et al., 2013;Lacoste et al., 2014;[START_REF] Whiteus | Perturbed neural activity disrupts cerebral angiogenesis during a postnatal critical period[END_REF]. However, the extent of these changes at the whole brain scale is unclear because the focus has so far been either on the capillary network in a small region (Harb et al., 2013) or on arterial network in a large region [START_REF] Lu | Netrin-1 hyperexpression in mouse brain promotes angiogenesis and long-term neurological recovery after transient focal ischemia[END_REF].

As a model of adult plasticity, we examined the nature of longterm changes to the vascular network 3 weeks after an ischemic stroke provoked by electrocauterization of the middle cerebral artery, a procedure that creates a focal ischemia in the dorso-lateral cortex. To avoid damaging potentially fragile vessels, we prepared drop-fixed brains from control and stroked animals, which we stained for podocalyxin, Acta2, and mouse immunoglobulins (as in Figures S5A andS5B). We first examined the middle cerebral artery (Figure 6A). The branch connecting the ventral circle of Willis to the dorsal arbor was missing (a stub of the branch is still visible ventrally), while the arterial dorsal arbor originating from the former degraded branch was preserved, hinting to the maintenance of the dorsal perfusion routes. When compared the vascular graphs of 3 control and 3 infarcted animals (Figures 6B and6C; Video S3). The center of ischemia in the supplemental somatosensory cortex (SSs, 2 brains) or in the nearby mouth somatosensory cortex (SSp-m, 1 brain) was noticeable from the scans due to a higher tissue autofluorescence and the presence of cellular debris. The capillary branch point densities outside of the ischemic area was not significantly different from the controls when averaged as voxelized densities (Figures 6B, 6C, and6E). However, a closer inspection of the center of ischemia and its surrounding volume revealed a striking reorganization of the capillaries so they were more oriented toward the ischemic center (Figure 6D). Orientation color-coding of the vessels from the graph highlighted the reduction in radial orientations (in green) in proportion to the vessels oriented along the cortical plane (in blue/red, Figure 6D). Following this observation, we quantified the vessel orientations between groups by voxelizing accumulated orientations. Accumulated orientations were different in the cortical regions of the stroked group around the ischemic center with a loss of the radial component (Figure 6F). The vessel orientation remodeling extended across several cortical regions, but didn't propagate to subcortical regions. This result suggests that density-preserving remodeling of vessel orientations is a strong factor in adult vascular plasticity in a model of induced ischemia.

The vascular topology can adapt to changes in neuronal activity levels early in life (Harb et al., 2013;Lacoste et al., 2014;[START_REF] Whiteus | Perturbed neural activity disrupts cerebral angiogenesis during a postnatal critical period[END_REF]. Long-range effects of a loss in specific sensory inputs on the organization of the vasculature in unrelated brain regions have not been studied yet. Because we previously noticed that the auditory system has higher densities of branching vessels than neighboring regions at most relays (Figures 4A and4B), we compared the vasculature of congenitally deaf mice to the one of control mice of the same background (FVB). Otoferlin (otof), defective in a recessive form of profound congenital deafness [START_REF] Roux | Otoferlin, defective in a human deafness form, is essential for exocytosis at the auditory ribbon synapse[END_REF][START_REF] Yasunaga | A mutation in OTOF, encoding otoferlin, a FER-1-like protein, causes DFNB9, a nonsyndromic form of deafness[END_REF] acts as a calcium sensor for synaptic vesicle release at the first relay synapse of the auditory pathways [START_REF] Michalski | Otoferlin acts as a Ca 2+ sensor for vesicle fusion and vesicle pool replenishment at auditory hair cell ribbon synapses[END_REF]. In Otof À/À mice, sound-evoked neuronal activity is absent but spontaneous activity during embryonic development is undisturbed, so that the initial neuronal projections still develop normally [START_REF] Mu ¨ller | Topographic map refinement and synaptic strengthening of a sound localization circuit require spontaneous peripheral activity[END_REF].

We generated graphs from 3 control and 3 Otof À/À mice 2 months of age. We did not observe a change in the number of penetrating arteries in the primary auditory cortex (Figure 7A), as previously reported [START_REF] Adams | The pial vasculature of the mouse develops according to a sensory-independent program[END_REF]. We compared then the branch densities of all vessels from averaged registered voxel maps obtained. There was a significant (p < 0.01), though small, reduction of branch point density in the Otof À/À auditory cortex as assessed from the voxel map. Strikingly, the p value map revealed an increase in branch densities in the somatosensory barrel cortex (Figure 7B), as well as in the visual cortex, suggesting a cross-modal compensation for the loss of auditory input. Intriguingly, a few integrative regions, such as the peri-and ectorhinal areas, strongly connected to sensory cortices also showed a significant increase in branch point densities in the deaf Otof À/À mice. Cell 180,12345678910111213141516February 20,2020 9 We next explored the subcortical relays and found that most of the nuclei of the central auditory pathways in Otof À/À mice had lower branching point densities than in control mice. These nuclei included the cochlear nucleus, superior olivary complex, nucleus of the lateral lemniscus and inferior colliculus (Figures 7C and7D). Although the barrel cortex in Otof À/À mice had much higher branch densities than in control mice, the somatosensory relays to the barrel cortex were not affected (Figure 7D), suggesting that cross-modal vascular plasticity was only present in the cortex. Finally, we examined whether the connectivity of capillaries was affected in Otof À/À mice by grouping vessels into sub-domains based on the strength of their connections (Figure 7E). We found that although the number of branches was significantly different between the groups in the layers 5 and 4 from the auditory and barrel cortices, respectively, the number of domains was unchanged (Figures 7F and7G). This finding suggests that the densification of the capillary network may not change vascular territories. Our findings extend the report from [START_REF] Whiteus | Perturbed neural activity disrupts cerebral angiogenesis during a postnatal critical period[END_REF] to other brain regions, showing the

DISCUSSION

We present a strategy to generate labeled, region-annotated graphs of the mouse brain vasculature. There has been a recent surge of interest in understanding how neuronal function interacts with the cerebral blood flow and how vascular dysfunctions potentiates some pathological states, such as neurodegenerative diseases [START_REF] Iadecola | The pathobiology of vascular dementia[END_REF][START_REF] Kisler | Cerebral blood flow regulation and neurovascular dysfunction in Alzheimer disease[END_REF]. Addressing these problems would clearly benefit from the development of image processing pipelines able to cope with the demands of large-scale 3D fluorescence microscopy. Apart from the particular antibody cocktail, our method is not specific to the brain vasculature. Also, other tubular structures could be mapped: nephrons, lung bronchioles, seminiferous tubules, lymphatic system, etc. Indeed, we purposely designed the pipeline on a complex set of principled filters instead of deep neural networks (reserving the latter just for the binary tube-filling application), facilitating its tuning to variations in acquisition parameters or other types of biological objects, adaptations that may be more difficult with a trained network (Belthangady and Royer, 2019).

We validated the combination of 3 markers (podocalyxin, Acta1, and vWF) to discriminate veins, arteries, and capillaries. However, we did not use the data obtained from quadruple channels imaging in our analysis due to important chromatic shifts in light sheet focus and curvature across such a wide range of wavelengths. Therefore, our method would benefit from next-generation light sheet microscopes that improve on resolution and/or acquisition speed [START_REF] Chakraborty | Light-sheet microscopy of cleared tissues with isotropic, subcellular resolution[END_REF][START_REF] Pende | High-resolution ultramicroscopy of the developing and adult nervous system in optically cleared Drosophila melanogaster[END_REF][START_REF] Tomer | Advanced CLARITY for rapid and high-resolution imaging of intact tissues[END_REF][START_REF] Voigt | The mesoSPIM initiative: open-source light-sheet microscopes for imaging cleared tissue[END_REF]. Because access to such microscopes is still limited for most laboratories, we propose here a framework for precise reconstructions with high-throughput and accuracy on commercially available systems. Our method fails to correctly segment and label the pial and meningeal vasculature at the surface of the brain, due to the close apposition of veins and arteries in this planar network. Increasing the imaging resolution may not be sufficient to address errors in arterio-venous junctions visible at the surface of our graphs. Combining our method with high resolution imaging of the blood flow with fast ultrasounds (fUS) [START_REF] Errico | Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging[END_REF][START_REF] Hingot | Microvascular flow dictates the compromise between spatial resolution and acquisition time in Ultrasound Localization Microscopy[END_REF] could solve this issue by merging flow-based reconstructions of large vessels with our graphs, while at the same time adding speed and direction information on vessel branches.

The throughput and reliability of the method in generating vascular graphs implies broad possibilities for expanding and revisiting the role of the vasculature in neural functions and pathology. The pipeline lowers the technical barrier to address difficult questions, such as the mechanisms controlling the plasticity of the vasculature in the postnatal brain. Our analysis of a model of ischemic stroke revealed a deep re-orientation of the capillaries toward the anoxic site, an effect that has been documented in other models, such as traumatic brain injury [START_REF] Cao | Three-dimensional imaging of microvasculature in the rat spinal cord following injury[END_REF][START_REF] Sosa | Low-level blast exposure disrupts gliovascular and neurovascular connections and induces a chronic vascular pathology in rat brain[END_REF]. However, although the addition of new endothelial cells to cerebral vessels has been reported in models of ischemic strokes [START_REF] Liu | Vascular remodeling after ischemic stroke: mechanisms and therapeutic potentials[END_REF][START_REF] Zhang | Bone marrow-derived endothelial progenitor cells participate in cerebral neovascularization after focal cerebral ischemia in the adult mouse[END_REF], we did not measure an increase in branch point density 3 weeks after the stroke away from the site suggesting that the reorganization of the vasculature following hypoxia at a distance from the infarct may involve density-preserving remodeling.

In contrast, we measured opposite effects in a model of activity-dependent developmental plasticity. In congenitally deaf mice, we noted a change in the vascular density that does not affect the size and numbers of capillary territories. Sensory deprivation is known to trigger in the cortex cross-modal compensation in neuronal activity levels, usually by lifting longrange inhibition [START_REF] He | Distinct sensory requirements for unimodal and cross-modal homeostatic synaptic plasticity[END_REF]. We could measure this effect by Fos mapping after whisker trimming, which led to an increase in Fos + cells density in the auditory cortex (Renier et al., 2016). Much is left to do to characterize the rules of activity-dependent vascular plasticity in the brain, and our tool should facilitate this line of inquiry, because it enables the detection of subtle effects over distant regions. Moreover, the plasticity of the vascular network is intimately linked with neuronal developmental and adult plasticity [START_REF] Paredes | Neurovascular Communication during CNS Development[END_REF]. Our tools will help illuminate topological remodeling of both networks in the developing and adult brains.

To further elucidate the plasticity of the vascular network topology, our method enables the molecular annotation of vessel branches on top of the reconstructed network. This annotation will be invaluable to quantify heterogeneity in marker expression in the endothelium and surrounding cells. For instance, a novel population of Pdgfra + fibroblasts located in the Virchow-Robin perivascular spaces of arteries and veins has recently been described [START_REF] Vanlandewijck | A molecular atlas of cell types and zonation in the brain vasculature[END_REF]. Interestingly, these cells express the marker LAMA1 when associated with the vasculature, which could provide an extra marker to discriminate veins in combination with Acta2. As an extension, and because it is based on immunolabeling with/without intracardiac fixative perfusion, our method could also be used to map the expression of markers of pericytes, inflammation, Aß deposits, hypoxia or pathogens on vessel branches, the location of resident immune cells, and microglia close to the vasculature throughout the brain.

Finally, the pipeline could also help investigate the relationship between neuronal networks dysfunctions and metabolism. The implication of a vascular component has been suggested for several disorders affecting higher cognitive functions such as autism [START_REF] Azmitia | Persistent Angiogenesis in the Autism Brain: An Immunocytochemical Study of Postmortem Cortex, Brainstem and Cerebellum[END_REF], schizophrenia [START_REF] Hanson | Theories of schizophrenia: a geneticinflammatory-vascular synthesis[END_REF], or obesity [START_REF] Dorrance | The effects of obesity on the cerebral vasculature[END_REF]. However, the technical demands and time-cost for evaluating the link between (E) Averaged voxelized branch densities measured from the registered graphs of controls or post-ischemic brains (3 weeks), and the p value map (0.01) of their differences (Mann-Whitney test on voxels, n = 3). No statistically significant changes are observed between the groups. (F) Voxelized branch orientations, 1 representative example of 3 shown. In the control brain, a mix of all orientations can be seen in the cortex, whereas in the postischemic brain, a decrease of radial orientations (green) is seen far from the ischemic site (dotted circle). Scale bars represent 1 mm, except in (D) (300 mm).

vascular development/remodeling and neuronal function is daunting. Our technique provides a new powerful tool of investigation by making large-scale comparisons of vascular topology in disease models accessible.

STAR+METHODS

Detailed methods are provided in the online version of this paper and include the following:

d KEY RESOURCES

STAR+METHODS KEY RESOURCES TABLE LEAD CONTACT AND MATERIALS AVAILABILITY

Further information and requests for resources and reagents should be directed, and will be fulfilled by the Lead Contact, Nicolas Renier (nicolas.renier@icm-institute.org). This study did not generate new unique reagents. EXPERIMENTAL MODEL AND SUBJECT DETAILS

Animals

All procedures followed the European legislation for animal experimentation (directive 2010/63/EU). Animal manipulations were approved by the institutional Ethical Committee (Project Ce5/2016/3996). Adult male and female mice were used in the study of C57Bl/N and FVB backgrounds, obtained from Janvier Labs. All the animals were in between 8 and 12 weeks old when the experimental procedures were initiated. Animals were housed in group, in ventilated racks with free access to food and water. Temperature, humidity and photoperiod were controlled.

METHOD DETAILS

Perfusion and tissue processing All animals were sacrificed by Pentobarbital overdose (200mg/100 g of animal). Intracardiac perfusion was then performed with a peristaltic pump (Gilson, USA). Blood was washed by infusing 30mL of cold PBS followed by tissue fixation with 30mL of cold 4% Paraformaldehyde (Electron Microscopy Sciences, USA) diluted in PBS. Then the brains were quickly dissected taking special care to preserve intact the structure. Finally, the brains were postfixed for 2 hours at room temperature by immersion in 4% Paraformaldehyde, and stored in PBS at 4 C until further processing.

Samples staining and iDISCO+ clearing

Whole brain vasculature staining was performed following the iDISCO+ protocol previously described (Renier et al., 2016) with minimal modifications. All the steps of the protocol were done at room temperature with gentle shaking unless otherwise specified. All the buffers were supplemented with 0,01% Sodium Azide (Sigma-Aldrich, Germany) to prevent bacterial and fungi growth. Perfused brains were dehydrated in an increasing series of methanol (Sigma-Aldrich, France) dilutions in water (washes of 1 hour in methanol 20%, 40%, 60%, 80% and 100%). An additional wash of 2 hours in methanol 100% was done to remove residual water. Once dehydrated, samples were incubated overnight in a solution containing a 66% dichloromethane (Sigma-Aldrich, Germany) in methanol, and then washed twice in methanol 100% (4 hours each wash). Samples were then bleached overnight at 4 Ci n methanol containing a 5% of hydrogen peroxide (Sigma-Aldrich). Rehydration was done by incubating the samples in methanol 60%, 40% and 20% (1 hour each wash). After methanol pretreatment, samples were washed in PBS twice 15 minutes, 1 hour in PBS containing a 0,2% of Triton X-100 (Sigma-Aldrich), and further permeabilized by a 24 hours incubation at 37 CinPermeabilization Solution, composed by 20% dimethyl sulfoxide (Sigma-Aldrich), 2,3% Glycine (Sigma-Aldrich, USA) in PBS-T. In order to start the immunostaining, samples were first blocked with 0,2% gelatin (Sigma-Aldrich) in PBS-T for 24 hours at 37 C. The same blocking buffer was used to prepare antibody solutions. A combination of primary antibodies targeting different components of the vessel's walls was used to achieve continuous immunostaining. Antibodies to Podocalyxin and CD31 were combined to stain the full capillary net and large veins, while an anti-Alpha Smooth Muscle Actin (aSMA or Acta2) antibody was used to label the artery's wall (antibodies' references and concentrations are provided in the Key Resources Table ). Primary antibodies were incubated for 10 days at 37 C with gentle shaking, then washed in PBS-T (twice 1 hour and then overnight), and finally newly incubated for 10 days with secondary antibodies. Secondary antibodies conjugated to Alexa 647 were used to detect Podocalyxin and CD31, while arteries were stained with secondary antibodies conjugated to Alexa 555. After immunostaining, the samples were washed in PBS-T (twice 1 hour and then overnight), dehydrated in a methanol/water increasing concentration series (20%, 40%, 60%, 80%, 100% one hour each and then methanol 100% overnight), followed by a wash in 66% dichloromethane -33% methanol for 3 hours. Methanol was washed out with two final washes in dichloromethane 100% (15 min each) and finally the samples were cleared and stored in dibenzyl ether (Sigma-Aldrich) until light sheet imaging.

Blood retention

The blood retention experiments were designed to retain as much blood in the vessels as possible while optimizing fixative diffusion in the tissue. These experiments involve a four-day process of dissection and post fixation steps. First, the mice were euthanized with CO2 and their skin was carefully dissected from the neck to the snout, avoiding injury to the surrounding blood vessels, especially the jugular veins. Following the dissection, the animals were placed in 4% Paraformaldehyde (PFA) for post-fixation, at 4 C, overnight. The next day, the interparietal bone was dissected, the snout was removed with surgical scissors, and the animals were again placed in 4% PFA for post-fixation, at 4 C, overnight. On the third day, the brain was dissected and placed again in 4% PFA for a final postfixation step, at 4 C, overnight. On the fourth day, the brains were washed 4-5 times with 1X phosphate-buffered saline (PBS) and kept at 4 C before further processing.

Ischemic stroke

Ischemic stroke was induced by middle cerebral artery (MCA) electrocoagulation as previously described [START_REF] Llovera | Modeling stroke in mice: permanent coagulation of the distal middle cerebral artery[END_REF][START_REF] Tamura | Focal cerebral ischaemia in the rat: 1. Description of technique and early neuropathological consequences following middle cerebral artery occlusion[END_REF]. Briefly, animals were deeply anesthetized with isoflurane 5% and maintained with 2.5% isoflurane in a 70%/30% mixture of NO2/O2. Mice were placed in a stereotaxic device, the skin between the right eye and the right ear was incised, and the temporal muscle was retracted. A small craniotomy was performed, the dura was excised, and the middle cerebral artery (MCA) was e2 Cell 180, 1-16.e1-e13, February 20, 2020 exposed. The MCA was electro-coagulated before its bifurcation with bipolar forceps, under irrigation, then the craniotomy, muscle and skin were closed and sutured.

Light sheet imaging

As a compromise between imaging speed and isotropic axial resolution, we imaged with a 4X 0.35NA objective cropped elongated field of view (600 3 2200mm) covering the narrow waist of the light sheet at 1.63mm/pixel of lateral resolution. Increasing the magnification didn't improve the axial resolution (Figure S3C). The planes spacing and cropping factor where optimized to avoid fusion of nearby vessels (Figures S3D-S3F). We observed that this tiling configuration led to a small anisotropy in the measured capillary diameters (Figures S3G and S3H). While shifts between the channels due to chromatic aberration were observed (Figure S3A), corrections were not necessary as the arterial binary image mask still overlap with the vessels binary image mask. This approach yields a near isotropic resolution, while striving to maintain an acquisition time under 1 day per brain (6h per hemisphere per channel) by expanding the width of the FOV over the ideal Rayleigh distance:

Z r = 2p$u 2 0 l
where u 0 = ðn $l =p NAÞ is the beam waist. At this NA and 647nm wavelength, Z r is 100mm. Tiling the cropped 4X FOV generates between 72 and 91 tiles per hemisphere, between 2900 and 3200 planes at 1.6mm spacing, which yields about 200Gb of uncompressed image files per channel per hemisphere. A reference channel for the registration to the annotated atlas using the sample autofluorescence was acquired at 5mm/pixel. As we are using the very limits of what enables a homogeneous image, it is extremely important that the microscope light sheet beam is well centered in all directions: y position, x,y tilt and planarity.

The acquisitions were done on a LaVision Ultramicroscope II equipped with infinity-corrected objectives. The microscope was installed on an active vibration filtration device, itself put on a marble compressed-air table. Imaging was done with the following filters: 595/40 for Alexa Fluor-555, and À680/30 for Alexa Fluor-647. The microscope was equipped with the following laser lines: OBIS-561nm 100mW, OBIS-639nm 70mW, and used the 2nd generation LaVision beam combiner. The images were acquired with an Andor CMOS sNEO camera. Main acquisitions were done with the LVMI-Fluor 4X/O.3 WD6 LaVision Biotec objective. The microscope was connected to a computer equipped with SSD drives to speed up the acquisition.

The brain was positioned in sagittal orientation, cortex side facing the light sheet, to maximize image quality and consistency. Indeed, while horizontal orientation is possible and enables whole-brain scanning in one acquisition, it is impossible here, as we use 1 light sheet illumination to increase the axial resolution to the maximum (combining 3 or 6 light sheets would be very detrimental to the resolution). Using 1 light sheet illumination in horizontal scan reduces illumination consistency because of 1) Excessive shadows from the white matter tracts 2) Light scattering at the midline of the tissue (see [START_REF] Ariel | UltraMicroscope II -A User Guide (University of North Carolina at Chapel Hill[END_REF] for a complete guide and explanations).

A field of view of 400 3 1300 pixels was cropped at the center of the camera sensory. The light sheet numerical aperture was set to the maximal NA (0.1). Beam width was set to the maximum. Only the center-left light sheet was used. Laser powers were set to 100% (639nm) or 40% (561nm). The center of the light sheet in x was carefully calibrated to the center of the field of view using capillaries, which should appear as points, not lines, when the light sheet is centered. z steps were set to 1.6mm. Tile overlaps were set to 10%. The acquisition routine was set to 1) Z-drive -> Save ome.tif stack 2) Filter change -> Z-drive -> Save ome.tif stack 3) Change X position -> repeat 1,2 12 times 4) Change y position -> repeat 1,2,3 6 times. The whole acquisition takes about 14h per hemisphere.

At the end of the acquisition, the objective is changed to a MI PLAN 1.1X/0.1 for the reference scan at 488nm excitation (tissue autofluorescence). The field of view is cropped to the size of the brain, and the z-steps are set to 6mm, and light sheet numerical aperture to 0.03 NA. It is important to crop the field of view to the size of the brain for subsequent alignment steps. After this acquisition (about 3 minutes), the brain is flipped to the other hemisphere and upside-down, so that the cortical side still faces the left light sheet.

The signal-to-noise (SN) ratio was calculated as:

SN = max A ðÞ À m B s B
where m B is the mean of the background signal, max(A) is the peak of signal in the weakest capillaries and s B is the standard deviation of the background signal. Full Width at Half-Max (FWHM) where calculated by fitting Gaussian functions to the profiles of capillaries taken in different location and orientations of the scans.

Computing resources

The data were automatically transferred every day from the acquisition computer to a Lustre server for storage. The processing with TubeMap was done on local workstations, either Dell Precision T7920 or HP Z840. Each workstation was equipped with 2 Intel Xeon Gold 6128 3.4G 6C/12T CPUs, 512Gb of 2666MHz DDR4 RAM, 4x1Tb NVMe Class 40 Solid State Drives in a RAID0 array (plus a Cell 180, 1-16.e1-e13, February 20, 2020 e3 separate system disk), and an NVIDIA Quadro P6000, 24Gb VRAM video card. The workstations were operated by Linux Ubuntu 18.04LTS. TubeMap was used on Anaconda Python 3.6 or Python 2.7 environments.

WobblyStitcher

Our method to acquire isotropic data via light-sheet imaging produces a large number of tiles that have to be aligned and stitched. For accurate reconstruction of the vasculature network at the vasculature capillary level a precise and robust alignment of large datasets is needed.

Existing stitching tools for large datasets could not generate accurate stitching results for our datasets due to misalignments arising from the placement methods (TeraStitcher) or large data reading and writing (BigStitcher). Even for well-placed tiles, rigid stitching produced object duplications in sub-regions along the z-axis indicating x-y movements of the images along this axis. Those movements likely originate from imprecise stage movements and were found on different microscopes and are prevalent in other datasets from other labs (personal communication from Pablo Ariel, University of North Carolina). To solve these problems, we developed WobblyStitcher, an open-source software tool for robust and non-rigid stitching of terabyte datasets.

WobblyStitcher uses three main steps to stitch the data (Figure 1H): (i) alignment, (ii) placement, and (iii) stitching. In the alignment step the non-rigid displacements between images are estimated along with additional quality measures of the alignment reliability. In the placement step a combination of tracking and optimization tools is used to robustly position individual tiles within a global image frame. Finally, the stitching step uses 3D interpolation methods in the overlapping regions to compose the final volumetric image.

The 3D non-rigid alignment is done in three sub-steps: (a) rigid alignment of the tiles along the axis orthogonal to the light sheet plane (z-axis), (b) placement of the tiles along this axis, followed by (c) non-rigid alignments within each plane (x-y-planes).

While a full non-rigid 3D alignment is possible by employing for example non-linear registration methods (e.g., elastix) [START_REF] Klein | elastix: a toolbox for intensity-based medical image registration[END_REF] we did not find this necessary for our datasets. Our non-rigid alignment is targeted toward light-sheet images and optimized for fast processing of terabyte-sized datasets. However, the modular architecture of WobblyStitcher is designed for extensions including full 3D non-rigid stitching and will be published elsewhere.

Rigid alignment

In the first step, neighboring tiles are rigidly aligned along the z-axis. To speed up this step, a max-intensity projection (MIP) is calculated along the axis that connects the two tiles and the alignment is performed on the resulting 2d images. While WobblySticher implements full 3d rigid alignment, we did not find this necessary for our datasets in this step and opted for the faster MIP version.

The n-dimensional rigid alignment is done via minimizing the mean square difference between the pixel intensities of two images I i and I j over all possible n-dimensional displacements s in the overlap region O ij ðsÞ as O ij ðsÞ is the number of pixels in the overlap. We set S ij ðsÞ = N if there is no overlap, i.e., N ij ðsÞ = 0. This estimate allows for a fast implementation via Fast Fourier Transform.

s à ij = argmin
WobblyStitcher implements correlation or normalized cross-correlation measures to estimate the shifts. In addition, other more sophisticated measures to estimate the displacements exists [START_REF] Klein | elastix: a toolbox for intensity-based medical image registration[END_REF] but are not implemented in WobblyStitcher yet, as we did not find those necessary to align our datasets.

Differences in high intensity voxels will contribute strongly to the error measure S above and can induce unwanted fluctuations, thus WobblyStitcher gives an option to preprocess the images before calculating the alignment. Arbitrary preprocessing routines can be passed. For the brain vasculature preprocessing includes clipping at specified intensity values above a threshold q h and below q l as well as normalizing the images by subtracting the mean and dividing by the standard deviation in the overlap regions.

In addition, the alignment estimates s à ij can be corrupted in regions without or very little foreground signal as background noise is aligned in this situation. WobblyStitcher provides the option to pass a validation routine to the alignment and to measure the quality of the alignment to prevent this problem. For the brain vasculature the validation is done on the raw data by requiring that the number of foreground pixels in the overlap region exceeds a certain fraction or minimal number of required valid pixels, i.e È q b <x<q t j x ˛Oij ð0Þ É Rn v

where q b and q t are intensity thresholds below or above pixels are considered as invalid for alignment.

To measure the quality of the alignments the following measure is used

q à ij = ( ÀS s à ij for valid alignments ÀN else:
e4 Cell 180, 1-16.e1-e13, February 20, 2020

Optimal rigid placement

Having measured the displacements between all tiles along the z-axis a placement along this axis is performed. For this we use globally optimal placement strategy to find the position p i for each tile i by minimizing the error term

E = X ði;jÞ˛P À p i + s ij À p j Á 2
with respect to the positions p i . Here P is the set of pairs of indices ði; jÞ of all neighboring tiles for which the alignment quality q à ij > q q exceeded a certain quality threshold q q . The resulting set of equations is solved for the positions via the Penrose pseudo inverse. The equations can be solved separately in each coordinate dimension to speed up computation of full 3d rigid placements or to calculate the optimal alignment along a single axis. More importantly, thresholding the alignment quality can force the alignment graph (nodes given by the tiles and edges indicating valid alignment pairs of tiles) to become disconnected. WobblyStitcher accounts for this by optimizing each connected component individually. For the brain vasculature the quality threshold is set to q q = À N in the z-alignment step only separating empty tiles from the brain sample.

Non-rigid alignment

In the final step every pair of tiles is aligned non-rigidly in the x-y planes via a series of computations:

First, the overlap region for each pair of neighboring and already z-aligned tiles is computed. The overlap region is validated via an optional validation routine passed to WobblySticher. For the brain vasculature we use foreground pixel counting as described above. The overlap regions are also preprocessed as a whole using normalization.

Second, for each z-value in the overlap region slices of the tiles in the x-y plane are taken and the alignment error S defined above computed for each slice. For the vasculature each slice is validated again separately using raw data slices and a higher threshold for the number of required valid pixels n v .

Third, the best displacements along the z-axis of the two tiles are determined via tracking. The reason for this step is to achieve robustness against errors that can occur as multiple local minima can arise in the error landscape S in datasets with repetitive structures near the overlap region. When considering global minima alone, this can create jumps in the displacements and errors in the alignments. WobblyStitcher first creates a list of all local minima in the error landscape S (excluding pixel at the overlap border) together with their quality measure q for each x-y-slice. Invalid marked slices will separate continuous segments of valid slices along z. For each continuous segment linear programming [START_REF] Jonker | A shortest augmenting path algorithm for dense and sparse linear assignment problems[END_REF] is used to track the local minima resulting in potential displacement paths through z. To find the full path in each segment the longest potential path with best quality measure is selected and all trajectories within that z-range of paths removed. Subsequently, the next longest path with minimal alignment error is selected etc. Valid slices that could not be assigned a displacement in this way will be marked as untraced and invalid.

Fourth, an optional smoothing is applied to the segments of valid displacements along the z-axis. For the vasculature we use convolution with a Bartlett window to smooth the segments.

The total result of this alignment procedure is that tiles become aligned in z and for each neighboring x-y-plane a robust displacement is determined together with markers for validity, quality and traceability.

Final Placement

In this second step the tiles are placed non-rigidly considering the validity and quality of the displacement measure.

WobblyStitcher first use the quality based global optimal alignment method described above to place the individual x-y planes in each z-slice of the entire dataset. For each slice this results in optimal displacements of the connected components. We will refer to each of these connected components in a single plane as a cluster for simplicity in the following.

In the second placement step, clusters are aligned optimally in the entire image. To achieve this the connectivity structure between the clusters is determined first by constructing a graph with clusters as nodes and edges between two clusters if they are from subsequent z-planes and overlap in the x-y-plane. The connected components of clusters in this graph are then aligned, by considering all displacements between the clusters and optimizing the error function:

E = X s X i˛Cs X j˛C s + 1 X k˛C s;i XC s + 1;j ÀÀ p s;k + d s;i Á À À p s + 1;k + d s + 1;j ÁÁ 2
with respect to the unknown displacements d s;i for cluster i in slice s. Here C s is the set of clusters in slice s belonging to the connected component under consideration, C s;i is the index set of all tiles of the cluster i in slice s, and p s;k is the already optimized valid position of the slice s of tile k. The optimization is again performed via the pseudo inverse to solve for the displacements d s;i . This results in optimal cluster positions and consequently in optimal positions of the individual image planes.

In a post-processing step, segments of untraceable or invalid slices of tiles that are neighboring valid slices in the z-axis are positioned using linear interpolation between the positions of the valid slices. Finally, the displacements along the z-axis are smoothed again via convolution with a Bartlett window.

To speed up computation in large datasets WobblyStitcher implements the option to perform the alignment and placement only on a subset or sub-grid of z-planes and fill in the skipped slices via interpolation.

Cell 180, 1-16.e1-e13, February 20, 2020 e5

Stitching

In this final step the tiles are combined into a single image. Stitching is done in each z-slice separately by using the optimized positions of the individual tile slices and their shapes to determine the various regions of overlaps together with the tile slices contributing to each of those regions. In addition, for each pixel in each tile slice the distance to the border of that tile slice is determined and passed as a weight for each overlap region to the actual stitching function. WobblyStitcher implements a set of stitching functions, such as maximum, minimum or mean projections as well as interpolations using the pixel weights. Also custom defined stitching functions can be passed. For the vasculature we use a weighted mean using the distance to the border as a weight.

ClearMap 2.0 toolbox ClearMap 2.0 is designed for advanced and fast image processing of large (Terabyte) 3D datasets obtained from tissue clearing. It is a complete redesign of ClearMap (Renier et al., 2016). The modular source management allows a unified handling of various data sources (image files, binary files, memory maps, shared memory arrays, numpy arrays, GPU arrays, graph formats) and fast parallel 3D image and graph processing as well as interactive visualization of large 3d images.

In

TubeMap pipeline

To construct graphs form the vasculature network we designed the processing pipeline TubeMap in ClearMap 2.0 detailed in the following.

Overview: a pseudo-deconvolution step removes shadowing artifacts of large vessels (Figure 2D); an adaptive thresholding captures signal inhomogeneities between and within large filled vessels; an equalization step homogenizes signal gradients across the sample and boosts weaker vessels and capillaries (Figure 2E); a modified tube filter amplifies smaller sized capillary structures with low signal to noise ratio (Figure 2F) Clipping, masking and first binarization (step 1) In a first step of the processing pipeline the stitched raw volumetric image is clipped above intensities that can be unambiguously assigned to foreground pixels. Because of 'stray' or 'blur' artifacts from bright vessels (Figure 2D) this upper intensity is chosen conservatively and will only include very large and bright vessels. The voxels clipped in this way contribute the first foreground voxels to the final binary image.

Voxels outside the brain sample show a distinct lower intensity value than even background voxels within the brain. Thus, in parallel to the clipping of high intensities, low intensities are clipped below the sample background and the clipped pixels are designated as background. Non-background pixels are used to define a mask for the brain sample which is used throughout the processing pipeline to restrict all calculations to the brain sample voxels. Light-sheet artifact correction (step 2) Light-sheet microscopy introduced 'stripe' and 'shadow' artifacts (Figure 2C). While methods exist to correct for those artifacts (e.g., via Fourier analysis, or stripe detection), they have a high computational demand making them difficult to apply to images of TB size.

We thus developed a method for fast and efficient light-sheet artifact correction. The method uses the fact that the stripe artifact occurs along a predefined axis in the images as the light-sheet always enters the sample from the same direction, and we used a single fixed light sheet illumination. Moreover, background voxels typically show similar intensities along the stripe axis over a certain length scale.

Thus, for each voxel i we estimate the light-sheet stripe artifact intensity l i by calculating a predefined percentile p l of the voxel intensities in a region centered around i and highly elongated along the stripe-artifact axis (Figure 2C). The length of this region (or structuring element) for the percentile filter along the stripe artifact axis is chosen to be of the scale on which the stripe artifact intensity changes, while the width and depth are chosen to be below the size of the stripe artifacts cross-section. Subtracting this estimate from the image leads to good corrections of the stripe artifact except at voxels that are part of longer vessels aligned with the artifact axis. To prevent this type of voxels to be removed, we also estimate the local background intensity b i as the percentile p b of the voxel intensities in a square shaped region centered around i. The size of this region is chosen larger than the largest vessel structures in the brain sample. This background estimate is compared to the light-sheet artifact estimate and voxel intensity v i corrected according

v i / v i À minðv i ; minðl i ; fb i ÞÞ
where f is a factor allowing to adjust the background estimate. The result of this correction is shown in Figure 2C.

Processing speed is limited by the local percentile filters. We use our fast 3d rank filter library to speed up the computation. On top, the elongated shape of structuring element for the light-sheet estimate allows further optimization by shifting the structuring elements for local histograms estimation in the direction of the artifact axis. For background estimation we use large structural elements which allow the use of sub-sampling and interpolation. Median filter (step 3) In the third step the light sheet corrected image is 3d median filtered with a small structuring element to smooth voxel intensities but preserve edges. We use our fast sliding histogram 3d rank filter library for this step.

Pseudo deconvolution and second binarization (step 4a)

The data showed 'blur' or 'stray' artifacts producing 'halos' of high intensity voxels around bright and large vessels (Figure 2D). Without correction of these artifacts, high thresholds had the tendency to remove weaker and smaller vessels and capillaries connecting to the larger ones, while low threshold values have the tendency to join larger neighboring vessels, even when local or adaptive thresholding was used. While this artifact could be corrected via an appropriate deconvolution step, such a step would add a large computational overhead for our very large datasets.

Cell 180, 1-16.e1-e13, February 20, 2020 e7

We thus designed a 'pseudo deconvolution' step that corrects for the 'blur artifact' with much less computational demands. High intensity voxels are identified via a threshold and 'blured' using a 3d Gaussian filter. The blurred result is then subtracted from the original image and the image rectified while the values of the high intensity voxels are preserved (cf. Figure 2D).

The resulting 'deconvolved' image is then thresholded at a lower intensity level than the one used to determine the high intensity voxels in order to capture the structure of the bright vessels and the resulting binary added to the final binarized image. Adaptive threshold and third binarization (step 5a) In order to capture smaller and less bright vessels the deconvolved image is further subjected to a local histogram-based adaptive threshold and the result added to the final binarized image.

We use our local histogram sampling framework to locally apply the Ridler-Calvard method [START_REF] Ridler | Picture thresholding using an iterative selection method[END_REF]) that determines a threshold by separating the voxels of the image into two groups such that the threshold is midway between the mean intensities of these groups. Equalization and fourth binarization (step 4b) While the previous stream of processing and binarization steps (4a,5a) is designed to mainly capture larger or brighter vessels, in parallel we designed a second stream (4b, 5b) for the binarization of weaker and smaller vessels and capillaries. While our clearing and staining protocol is highly optimized to homogeneously label the vasculature, regional variations in the overall luminance exists and particularly smaller and weaker vessels in the deeper regions are often not detected by the first stream of binarizations.

In this step we thus apply a custom designed equalization filter to the median filtered image form step 3. The filter calculates for each voxel i a lower and upper intensity (l i and u i ) via a lower and upper percentile (p l and p u ) of the voxel intensities in a rectangular region R i centered around i. The region is chosen to be larger than the large vessel structures and our efficient histogram sampling framework is used to speed up computation. The voxel intensity is then normalized via

v i /f i v i with normalization factor f i = v i =l i u i =l i %m m=u i else
and m a maximal intensity value for the upper percentile. As a result the image shows a more homogeneous intensity distribution across the sample (Figure 2E). A fixed conservative threshold is applied to the normalized image and added to the final binary. Tube filtering and fifth binarization (step 5b) In a final step, a generalized Frangi tube filter [START_REF] Frangi | In Multiscale Vessel Enhancement Filtering[END_REF][START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF] to enhance smaller vessels and capillaries is applied. The filter first smoothes the image with a 3d Gaussian of a given scale and then calculates the sorted eigenvalues l 1;i R l 2;i R l 3;i of the Hessian matrix at each voxel i. The tubeness measure t i is then calculated according to

t i = 8 > > > > > > < > > > > > > :
l 3;i l 2;i l 3;i g 23 1 + l 1;i l 2;i ! g 12 l 1;i %0; l 3;i %l 2;i < 0 l 3;i l 2;i l 3;i g 23 1 À a l 1;i l 2;i ! g 12 a l 1;i l 2;i < 1; l 3;i %l 2;i < 0 0 else

(1) with parameters g 12 = 0:5,g 23 = 0:5 and a = 0:25 chosen to enhance tube like vessels but also account for their bending [START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF].

A threshold is applied to the tube filtered image to contribute the last part to the final binary image.

Binary filling

In this step, the combined binary image is subjected to a 3d binary filling operation. While the final image is binary and thus already smaller in size it is still a large array. Splitting the filling operation into subsets of the data would entail complex joining operations. TubeMap thus implements a parallel binary filling code based on flood filling from the border that operates on memory maps, allowing the binary filling of arbitrarily sized images not limited by memory.

Deep vessel filling

As antibodies only target vessels walls receptor proteins, larger vessels appear as empty tubes. Given that our method to extract the centerline that relies on topology-preserving 3D thinning, it is necessary to generate solid tubes on the binary mask. Therefore, we designed a method to detect and fill empty tubes to allow a correct centerline extraction via erosion. Vessel filling is a complex task on our data, as the shape, size and continuity of the vessel walls are variable. For this reason, we decided to use a deep convolutional neural network to solve this task.

e8 Cell 180, 1-16.e1-e13, February 20, 2020

Architecture

We based the initial network architecture on the DeepVesselNet architecture from [START_REF] Tetteh | Deepvesselnet: Vessel segmentation, centerline prediction, and bifurcation detection in 3-d angiographic volumes[END_REF], which detects the centerline from 3D scans of filled vessels. We iterated modification of this architecture to obtain the following: 2 maxpooling layers followed by 2 3D convolutional layers with dropout with kernel size 7 and 5 and 16 and 32 channels respectively. Then we added a depthwise separable convolutional layer of size 32 with kernel size of 3. Next, two new other convolutional layer are added followed by upsampling layers to get back to the initial input data size before the maxpool layers, and a last convolutional layer of size 2 followed by a center shifted sigmoid function to map output values between 0 and 1, and a softmax layer (Figure S5E).

The depthwise separable convolutional layer (desep conv) consists in a depthwise convolution followed by pointwise convolution. The idea is to increase the network performance while keeping the same amount of parameters.

Optimization function

Our loss function was composed of a classic binary cross entropy Loss to minimize the difference between the output and the ground truth fed to the network:

L = À BCE + CBL BCE = 1 N X N i = 1 y i log b y i + 1 À y i ðÞ log 1 À b y i ðÞ ð
where N is the number of samples and of a class balancing loss based on [START_REF] Tetteh | Deepvesselnet: Vessel segmentation, centerline prediction, and bifurcation detection in 3-d angiographic volumes[END_REF]:

CBL = L 1 + L 2 L 1 = À 1 jY + j X j˛Y + log Pðy i = 1ÞjX; WÞÀ 1 jY À j X j˛YÀ logðPðy i = 0Þ j X; W ! L 2 = À g 1 jY + j X j˛Y f + log Py j = 0 ðÞ X; W ÀÁ À g 2 jY À j X j˛Y f À log Py j = 1 ðÞ X; W ÀÁ g 1 = 0:5 + 1 jY f + j X j˛Y f + Py j = 0 X; W ÀÁ À 0:5 g 2 = 0:5 + 1 jY fÀ j X j˛Y fÀ jPy i = 1 j X; W ðÞ À 0:5 j where Y + , Y À , Y f + , Y fÀ
respectively represent the set of positive, negative, false positive and false negative labels.

Vessels account for a minority of voxels (about 15% of the data in our case). This makes the training memory and time consuming as large amounts of data are needed to reach convergence. Therefore, we used a loss which favorizes false positives and strongly penalizes false negative during the training, preventing the network to favor segmenting pixels as background.

Training

To generate a training vascular brain graph for training, we generated datasets comprised only of filled tubes. To obtain such datasets, we designed a preparation of whole head fixation, in order to retain the blood in the tissue (Figure S5A), and see the section ''blood retention'' for the preparation of these brains. After dissection, we immunostained the brain for circulating immunoglobulins [START_REF] Liebmann | Three-Dimensional Study of Alzheimer's Disease Hallmarks Using the iDISCO Clearing Method[END_REF], and complemented it with an immunostaining against podocalyxin to insure a better continuity of the labeled capillaries (Figure S5B), and we double-stained for Smooth Muscle Actin as before. We cleared and imaged the labeled brains with identical conditions to the endothelial wall-stained samples. Vascular graphs were generated using the same multi-path binarization pipeline (Figure 2A) and centerline extraction as all other datasets (Figure S5C). These graphs contain a natural distribution of vessel diameters, which we used to build training cubes of hollow tubes, for which the ground truth is the filled version (Figure S5D). We designed and trained the DNN on these data to fill hollow tubes into solid tubes.

We trained our network on 100*100*100 pixel wide blocks of artificially generated hollow vessels data. From the skeleton, we generated circular empty tubes as input training data and the filled counterpart as ground truth. To diversify the data, we also generated cubes from both the IgG/Podocalyxin dense capillary channel and from the Acta2 channel alone which has sparse arterial tubes. To diversify the vessel radii, we filtered out capillaries (which are already filled at this resolution) from the graph data to only keep vertices corresponding to vessels with larger radii. We also artificially increased the vessel radii on the sparse arterial graph to mimic underrepresented large vessels as they only account for a minority of the vessels in the dataset.

A preprocessing step was added, consisting in adding noise by setting randomly voxel values as foreground to mimic background dots that appears occasionally on binarized data, and train the network to ignore them.

We used a batch size of 8 and trained for 100 epochs.

Cell 180, 1-16.e1-e13, February 20, 2020 e9

Application to the vasculature data We used the network on the CD31 + Podocalyxin stained binarized scans (hollow vessels). The binary masks of 1 channel from 1 hemisphere represents approximately 100Gb and cannot be loaded to the network directly. Therefore we divided the binary masks of the generic vessels (Pdclx + CD31) channel in blocks of 500*500*500 voxel with a 100 voxel overlap in every direction.

We also filled the binary masks of the arterial channel. As Acta2+ vessels have larger radii and are sparse, we used blocs of 900*900*900 voxels with a 200 voxel overlap, and down sampled the blocks by a factor of 4.

Binary smoothing

In the final step we smooth the binary data, as rough surfaces of vessels can lead to artifacts downstream in the skeletonization of the vasculature network. We therefore developed a discrete topology based binary smoothing algorithm based on Ne ´meth et al. (2010).

For each voxel ithe local discrete topology in a 3x3x3 cube C i centered around i is considered and depending on the configuration the center voxel is updated as follows:

If the center voxel is foreground the center voxel is set to background if (i) there are less than three foreground voxels in the cube, (ii) the local topology matches the configuration C 1 or any of its rotations toward each of the 6 faces of the cube (iii) the topology matches C 2 or C 3 or any of the 4 rotations around the z-axis or any of those further rotated toward each face (i.e., 24 rotations), (iv) if the topology matches C 4 or any of its 12 rotations that rotate a fixed edge onto another one (v) the topology matches C 5 or any of its 8 rotations that rotate a corner onto another one.

If the center voxel is background the center voxel is set to foreground if either (i) the number of background pixels in the cube is less than 6, (ii) the inverse topology in which foreground is exchanged with background fulfills either of the conditions for the foreground pixel smoothing above, (iii) the center voxel has at least 3 neighbors in its 6-connected neighborhood and two of those voxels lie along one of the main axes.

The processing of the smoothing is optimized via assigning a unique topological index t to each possible topological configuration of the 3x3x3 cube around the center voxel and pre-calculating a lookup-table from this index to the resulting center voxel value according to the smoothing rules above. The index t for each voxel in the image is obtained by considering the binary cube configuration as a 27 bit representation of that index, i.e., by convolving the local cube with a kernel K. To further speed up the computation, we make use of the fact that the kernel K is separable and convolution with it can be calculated via three subsequent one-dimensional convolutions along the three axes.

Skeletonization and graph construction

In a final step of the TubeMap image processing pipeline the binarized image is converted to a graph. The graph construction is done in three steps: (i) the binary image is skeletonized, (ii) loose ends in the skeleton are detected and possible continuations to other loose ends detected via tracing in the non-binarized image data, (iii) the final skeleton is converted into a graph.

Skeletonization

In the first step of the graph construction the binary image is skeletonized. While a large number of skeletonization algorithms exist, they were not suitable for the size of our datasets with running times of over 7 days on our workstation. We thus implemented a fast skeletonization algorithm for terabyte volumetric images that runs in 45min. As skeletonization via thinning is context dependent the processing cannot be easily split into blocks of data. Instead, our algorithm is designed to handle TB sized data as whole. To increase the processing speed, we pre-calculate the thinning actions for all possible local topological configurations into a look-up table and use linear indexing of the image arrays. The algorithm uses 3d parallel thinning with 12 sub-iterations as described in (Pala ´gyi and Kuba, 1999) with a series of optimizations. In short, for each voxel i the algorithm decides if that voxel can be thinned away by considering the topology in the local 3x3x3 cube C i centered on i and using the rules defined in (Pala ´gyi and Kuba, 1999) (T1-T14). This process is done in parallel. As parallel thinning of surface voxels may lead to the disconnection of center lines or other topological changes in the final skeleton, the thinning is done in 12 sub-iterations 'attacking' the surface pixels from 12 different directions corresponding to the mid-points of the 12 edges of a cube.

To speed up processing on large images we implemented a series of additional optimizations. First, as the thinning is operating on foreground pixels only and their number significantly decays during the thinning iterations, calculation of the local topological configuration around each voxel is restricted to the foreground pixels only. In each thinning step only voxels that are 6-connected to the background are candidates for removal, so the calculation is further restricted to those border voxels. Second, the topological characterization around each voxel is done using topological indexing together with a pre-calculated look-up table that encodes the thinning rules as described above in the binary smoothing step. Third, all calculations are done on linear arrays using linear indexing, instead of 3d volumetric images. This speeds up addressing foreground pixels and reduces memory requirements by a factor of three. Fourth, all calculations, including candidate voxel detection and topological indexing are parallelized.

The skeleton is post-processed by removing foreground voxels of center lines with at least one endpoint (a foreground voxel with only one neighbor in the 27-neighborhood) and a length to a branch point (a foreground voxel with more than two neighbors) or another endpoint below a critical length. End-points and branch-points detection is efficiently done in parallel by calculating the number of neighbors via convolution with a 3x3x3 cube of ones around each foreground pixel.

e10 Cell 180,February 20,2020 End-point tracing In this optional step, the skeleton is further post-processed by trying to connect loose ends of the skeleton center-lines using tracing in the non-binarized tube-filtered image. To achieve this, the endpoints of the center-lines are detected as described above. For each endpoint an A*-search in the voxel space is performed to find a path with the lowest cost toward the closest endpoint. The cost C for a path P is defined as

C = X ði;jÞ˛P dði; jÞ t j
where ði; jÞ is a pair of neighboring voxels (27-neighborhood) in the path, dði; jÞ the Euclidian distance between the voxels, and t i is the tubeness measure defined in (1). If the cost for the path is below a threshold the path is added to the binary image. If the cost is too high, a path to the binary mask of the vasculature is traced using the same cost function. If the cost is below a threshold the path is added to the binary mask. The resulting binary is re-skeletonized as described in the previous step.

The purpose of this tracing step is to detect small vessels of weak intensity that were partly missed in the binarization process and thereby to reduce the number of open ends in the vasculature graph. For our datasets and with the above described equalization methods this step was usually not necessary (cf. Figure 2).

Graph construction

To facilitate the analysis of the vasculature network, the skeletonized binary image is turned into a graph representation that captures the topology of the vasculature network. A graph G = fV; Eg is a collection of vertices v i ˛V and edges e i;j ˛E between two vertices v i and v j . The vasculature network carries additional geometric information (e.g., vessel shape or radius) and can be accompanied by additional local data (e.g., expression levels of molecular markers). We thus consider graphs in which the vertices and edges hold a set of additional properties (p k ðv i Þ, p l ðe i;j Þ), such as spatial, geometric or molecular information.

Graphs extracted from the vasculature consist of millions to hundreds of millions of vertices. To enable the analysis of those large graphs, TubeMap provides a high-performance graph module based on the graph-tool library (Peixoto, 2014) and boost graph libraries https://www.boost.org/doc/libs/1_66_0/libs/graph/doc/index.html). In addition to a large number of graph manipulation and analysis routines, our module also provides graph classes that handle the spatial geometry and other annotational information, as well as the visualization of the graphs in 3d space.

In TubeMap, graph construction from the skeletonized image is done as follows: first, a raw graph is constructed by turning each foreground voxel of the skeleton into a vertex that also carries the positional information of the voxel. The vertices are then connected via undirected edges if they are neighbors (27-neighborhood) in the skeleton image. This process can result in local all-to-all connected cliques with more than 2 vertices at branch points of the skeleton. Thus, those cliques are identified and replaced by a single vertex with a coordinate position that is the mean of the clique vertices. Isolated vertices or small components not connected to the giant component are typically removed from the graph for downstream analysis.

In a second step, the raw graph is reduced to a branch graph consisting of vertices that are endpoints or branch-points only (i.e., vertices v i with edge degree dðv i Þs2). The vertices of this graph are then connected by edges if there is a path of vertices of degree 2 between them. The connection is done via tracing in the raw graph. In this process, the positional, geometric and other information attached to the vertices along a path between two branch-vertices is collected and attached as a property to the new edge connecting the two branch-points.

Vasculature graph geometry extraction

To capture the geometry of the vasculature network, TubeMap, besides storing the positional information about the center lines of the vessels also detects their shape. In particular, the radius of the vessel at each vertex or each point in a branch is estimated by taking the pre-processed equalized image and measuring the distance to the nearest voxel in which the intensity decayed by half. To avoid measuring radii within hollow tubes, the nearest voxel search is started outside the final binarized image in which hollow tubes are filled.

Vasculature vertex and branch labeling

TubeMap enables measuring the expression of other markers obtained via multi-color stainings and imagining along the vasculature graph. A set of tools is provided to measure those, including methods to measure the histogram, mean, maximum or minimum expression levels within a structural element centered on each voxel of the extracted graph or within a region described by the local extracted vessel geometry.

To measure the expression of the smooth muscle marker Acta2, we measure maximum expression levels around each point of the vasculature graph within a spherical region with a radius equal to the estimated radius of the vessel at that point. Branches are labeled positive for Acta2 if the majority of points in that branch have expression levels above a certain threshold.

TubeMap provides further routines to post-process edge and vertex labels based on geometric or other measures and label continuity. In particular, binary labels can be subjected to morphological operations acting on the graph topology. This includes binary morphological closing that can be used to fill gaps in a sequence of branch labels. In addition, labels of vertices or edges can also be traced along the graph topology according to arbitrary passed rules that act on the geometry or other vertex or edge properties, providing a generalized hysteresis thresholding operation for a given labeling. For the vasculature graph we use a one-step morphological closing as well as tracing based on the radial measure to post-process the Acta2 branch label. Reconstruction of Cell 180, 1-16.e1-e13, February 20, 2020 e11 the vasculature network was corrupted by dense vessel crossings as well as preparation artifacts on the brains cortical surface. Thus, to avoid errors in the artery and vein labeling the tracing of arteries and veins was stopped when the distance to the surface fall below a certain threshold. The distance of each branch to the brain surface in turn was calculated by an Euclidian distance transform from the background into the brain using the 3d brain atlas annotation image.

Vasculature graph atlas annotation

TubeMap provides functions to annotate data with labels from reference atlases, such as the Allen Brain Atlas. Annotation is done by acquiring an additional reference image in parallel to the other data channels and align this to the corresponding reference image of the atlas. We used an autofluorescence template from the ABA at 25mm generated from serial 2-photon tomography. To achieve the alignment, TubeMap first resamples the auto-fluorescent image of the sample to the resolution of the atlas reference and then aligns the image in 3d non-linearly. TubeMap integrates the elastix package functionality [START_REF] Klein | elastix: a toolbox for intensity-based medical image registration[END_REF] for the non-linear alignment. For the vasculature data a hierarchical estimation of b-spline transformations between the two reference images is used together with a cross-entropy measure to quantify the local alignment quality. For the vasculature data, the reference image is acquired in tiles directly after each tile of the other channels and thus the same alignment and stitching layout can be used to assemble the full reference image. To correct for small misalignments between reference and data images a rigid transformation between both is estimated using also a cross-entropy measure. The resulting joint transformation from data to sample reference to atlas reference is used to transform the coordinates of the graph points onto the reference frame of the atlas. In the final step, the transformed positions of the graph are used to extract the atlas annotation and any other atlas information.

The TubeMap graph module implements routines to extract statistics, sub-graphs or other information based on these annotations. The atlas annotation can also be used for 3d rendering of the vasculature graph.

QUANTIFICATION AND STATISTICAL ANALYSIS

Statistical analysis

All statistics were conducted with the Scipy Python library (https://scipy.org). Throughout the manuscript, data are represented as means with standard deviation. For comparisons of the means between two groups, non-parametric Mann-Whitney test were used (Figures 6E, 7A, 7B, 7D, and 7I). All sample sizes, precision measures and p values are indicated in the figure legends. In all instances, n represents the number of mice used. Further details on specific quantifications are given in the following sections.

Branch point density, branch density, branch length and tortuosity

Branch point density and branch density were defined as the number of branch points or branches per volume in each brain region, respectively. The volume of each brain region was estimated by summing the volume of all voxels belonging to that region using the reference atlas annotation.

The length of each branch was measured as the sum of the Euclidean distances between consecutive pairs of coordinates that define the branch as extracted from the raw graph.

Branch tortuosity was defined as the ratio between the length of the branch divided by the Euclidean distance between its endpoints.

Shortest paths to arteries and veins

Shortest paths from all branch points in the brain to the arteries and veins were computed via a graph-topology based distance transform algorithm. The graph-topology based distance transform was implemented in TubeMap via iteratively propagating artery or vein labels on the graph topology and at the same time keeping track of the iteration number at each propagation step (i.e., the shortest distance to the newly reached branch points). To consider the different lengths of each branch, the distance transform was computed on the raw graph instead of branch graph in which each branch is expanded to a chain of vertices corresponding to the branch length. While this introduces a small error when compared to the actual Euclidian lengths of the branches, it allows computing the shortest distances of all branch points at once, significantly speeding up the calculation and enabling a brain wide analysis. As a result, each branch point is assigned a distance to the closest artery branch as well as the closest vein branch. In a final step, only the distances at each original (non-expanded) branch point are kept. The shortest path length from arteries to veins running through a specified branch point was then calculated as the sum of those distances at the branch point under consideration.

Arterial domains

The capillary vasculature network is highly recurrent and redundant making it in principle difficult to identify domains that are 'fed' by a single arterial vessel. However, this analysis aimed at identifying those vasculature domains supported by arterial vessels based on shortest distances. In order to do so, we first identified isolated arterial vessels by labeling connected components in the arterial network (only the branches identified as arteries). To separate penetrating arteries into the cortex, the arterial network on the brain surface was removed before detecting the components. Each component was then given a unique label and all labels were iteratively propagated in parallel to the unlabeled capillary network. As a result, each capillary was assigned to its closest arterial component, defining the arterial domains.

e12 Cell 180, 1-16.e1-e13, February 20, 2020

Distances of vasculature graph features between cortical regions

The goal of this analysis was to analyze differences in the vascular organization and its layered structure between cortical regions. In order to do so, we used the ABA annotation to extract a set of features of the vasculature graph for each layer of each cortical region and compared those. Wee extracted the following five features: number of branching points (all vessels), number of arterial branching points, number of arterial end points, number of radially oriented vessels, and number of planar oriented vessels. The number of radially and planar oriented vessels was determined by first calculating the normalized orientation vector of the vessel. A vessel was counted planar (radial) oriented if the length of the orientation vector projected onto the planar plane (onto the radial axis) was larger than 0.7.

To finally compare cortical regions, we first normalized each feature resulting in a distribution of each feature across the 6 cortical layers. We then calculated the earth-movers distance (or ''Wasserstein distance'') between those distributions to obtain the overall distance between each pair of cortical layers. The earth-movers distance can be understood via the following picture: let each distribution be represented by a pile of sand. The distance between those distributions then simply is the work or cost to move the least amount of sand to transform one pile to another. Thus, the more the shapes of the distributions differ, the costlier it is to transform them, and the larger the distance between them.

Finally, after having obtained the distance matrix between all cortical regions, Ward clustering was performed to identify groups of cortical regions sharing similar distances to other regions.

Stochastic block model analysis

While the arterial domain analysis above was targeted at identifying sub-regions defined by arterial supply, in this analysis we aimed at identifying sub-networks or modules in the vasculature network based on the graph-topology alone. In order to do so we used an inference-based approach.

More precisely, we assumed that the adjacency matrix of the vasculature graph G is captured by a stochastic block matrix composed of B blocks with a number n s of nodes in each block s and a number e sr of random edges between blocks s and r with s; r˛f1; .; Bg (Holland et al., 1983). Probabilistic inference is then done on this model given the observed vasculature graph using the implementation described in (Peixoto, 2014). In short, as edges are assumed to be random between blocks, all graphs with the block structure C = ðfn s g; fe sr gÞ are equally probable. Thus the likelihood of observing a graph G that already has the block structure C is PðG j CÞ = UðCÞ À1 where UðCÞ is the total number of different graphs with block structure C. For a fixed block number B this likelihood can be maximized using a Markov Chain Monte Carlo algorithm that iteratively changes block memberships of individual nodes (Peixoto, 2014). In order to find the optimal block number B a minimum description length criterion is used to compare the different models (Peixoto, 2014[START_REF] Peixoto | Nonparametric Bayesian inference of the microcanonical stochastic block model[END_REF].

DATA AND CODE AVAILABILITY

Example datasets, links for ClearMap, ClearMap 2.0 and TubeMap codes download, tutorials and bench protocols are available online at https://idisco.info. Data are available on OSF: https://osf.io/sa3x8/?view_only=4427a838cbd0468c9fbad9cab465d866. The codes are also accessible from GitHub at https://github.com/ChristophKirst.

Cell 180, 1-16.e1-e13, February 20, 2020 e13 4 A-B Upper segments of the spinal cord (cervical to upper thoracic levels), labeled for CD31, Podocalyxin (blue) and acta2 (red). The anterior spinal artery is strongly positive for acta2, while the posterior spinal vein is also visible with a weaker acta2 signal. Collateralizations of the posterior spinal arteries are visible on each sides. Coronal projections of the reconstructed graphs highlight the arterial domains from the anterior spinal artery in the ventral horn and the posterior spinal arteries in the dorsal horns. Ingressions from the anterior and posterior veins are also annotated and visible in the midline projections. (C) Heatmap of the branch point density, example of a brachial level: higher densities of branch points are visible in the gray matter and pial surface compared to the white matter (D) Diffusion from arteries in the upper spinal cord. The left panel shows the distribution of the number of vessels across arterial distances (given as branch order from the closest artery). Most vessels are within 2 to 5 branches from the closest artery. This distribution is similar than the one calculated in the brain. The right panel shows the distribution of radii for each step away from an artery, which is also similar than in the brain, and shows larger radii (from venous components) from order 6 through 8.

IV.2. Summary of Article 2: The postnatal maturation and maintenance of the cerebral vasculature

The postnatal maturation and maintenance of the cerebral vasculature -draft: in writing 2023

From the previous study, it appears that vasculature is very heterogeneous across the brain regions and plastic to neuronal activity disruption in terms of orientation and density [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF].

We want to know what drives those topological heterogeneities, as it would give leads for developmental vascular patterning and remodeling mechanisms.

Despite the known existence of neurovascular coupling, how and whether or not they control each other's topology and function remains unclear, since there is no systematic correlation between neuronal and vascular density. Moreover, the adult vascular plasticity could never be properly investigated at a large scale because whole organ data acquisition has long been unachieved.

To i n v e s t i g a t e t h e s e , w e b u i l t a d e v e l o p m e n t a l a t l a s o f t h e g r o w i n g m o u s e b r a i n a t P 1 , P 3 , P 5 , P 7 , P 9 , P 1 2 , P14, and P21. We quantified the deformations induced by clearing by comparing MRI scans with light sheet autofluorescence. We then monitored the cerebral vascular network, total length, density, number of vessels, orientation distribution, across development. This enabled us to characterize two major developmental dynamics and the classification of the brain metaregions according to this dynamic, and modeling of vascular growth through different driving processes to understand how the cerebral vasculature caters and takes shape during the development and maintenance of its network.

Contributions: vessel filling DNN improvement, graph correction, improvement of the artery-vein labeling, development of analysis algorithms, characterization of the developmental waves, building the developmental atlas, alignment, and quantification of the structural differences between age-matched autofluorescent and MRI samples.

INTRODUCTION

The vasculature spreads as a densely intricate network through the brain to provide its neurons with oxygen and nutrients. Because of its low-capacity storage, the brain is very sensitive to any misadjustment between neuronal activity and tissue perfusion. Local changes in the blood flow are known to affect the neuronal patterns and conversely, neuronal activity locally modulates hemodynamics. ( [START_REF] Girouard | Neurovascular coupling in the normal brain and in hypertension, stroke, and Alzheimer disease[END_REF] [START_REF] Zlokovic | The blood-brain barrier in health and chronic neurodegenerative disorders[END_REF]) [START_REF] Ances | Coupling of changes in cerebral blood flow with neural activity: what must initially dip must come back up[END_REF] [START_REF] Iadecola | The pathobiology of vascular dementia[END_REF][START_REF] Kisler | Cerebral blood flow regulation and neurovascular dysfunction in Alzheimer disease[END_REF][START_REF] Zlokovic | Neurovascular pathways to neurodegeneration in Alzheimer's disease and other disorders[END_REF]. In [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF] it has been shown that, unlike a widespread belief, the vasculature is far from being homogeneous from one region to the other, suggesting that something drives the vascular organization regionally to ensure an accurate match between tissue perfusion capacity and local metabolic needs. Such organization strategy has been shown by [START_REF] Kennel | From whole-organ imaging to in-silico blood flow modeling: A new multi-scale network analysis for revisiting tissue functional anatomy[END_REF] in mice fat pad tissue. However, the neurovascular mechanisms are not fully characterized and very little is known about how vascular density and neuronal activity correlate. In (Blinder et al 2013), 2 photon was used to reconstruct the cortical angiome in the barrel fields, which is a highly heterogeneous region altering densely populated with neuron areas and axonal tracks. As no matching between vascular and neuronal patterns is observed, it is claimed that there is no link between vascular organization and neuronal activity. However, some other studies suggest the opposite to a certain extent. In (Lacoste et al 2014) show that early neuronal activity loss leads to hypervascularization, whereas [START_REF] Whiteus | Perturbed neural activity disrupts cerebral angiogenesis during a postnatal critical period[END_REF] show that early hyperactivity also leads to the same phenotype, suggesting, that development is a critical period where the vascular network's final organization is affected by local changes in the neuronal activity. Accurate analysis of whole organs often requires a reliable atlas of the corresponding structures. Therefore it becomes necessary to develop tools to study vasculature during the early days. However this is a very challenging task for the following reasons: from P1 to adulthood, the brain size spans across almost 1 order of magnitude and undergoes large-scale changes even in a few hours' time window. The different brain regions do not develop evenly, thus, to capture the exact timepoint where local changes happen many timepoints have to be analyzed. (Thompson et al 2014) and (Young et al 2021) propose atlases of the developing mouse brain, but only provide 2 postnatal timepoints: P4 and P14, acquired under MRI 2D acquisition, leading to incomplete and low resolved materials. To address this issue, we propose a pipeline to generate 3D atlas of the developing mouse brain at P1, P3, P5, P7, P9, P12, P14, and P21, an online resource to explore those data, and we characterize the evolution of brain vasculature at each timepoints.

RESULTS

Atlases generation pipeline

We acquired the 5×5×6 um/voxel light sheet auto-fluorescence from 10 brains at all the different timepoints.

Each hemisphere was upsampled to 25×25×25 um/voxel and then independently aligned using affine and bspline elastic transformation, resulting in an averaged image from 20 samples. The mouse brain not being considered lateralized, the hemispheres were mirrored and stitched together in order to get a full brain template. ITK-SNAP software was then used to manually annotate the template. 51 regions were annotated based on the Allen brain atlases. To eliminate the jerkiness induced by plan by plan drawing a laplacian smoothing was applied to the external mesh of each region (Fig. 1).

Data MRI vs iDISCO comparison

The iDISCO protocol used to generate the autofluorescence is known to induce anisotropic shrinkage of the initial tissue based on its initial densities and composition. To document those deformations, we compare the final volume of the different brain meta-regions in MRI and autofluorescence. Up to P9, little differences in volume have been witnessed. However, it shows that deformation due to skull pressure and CSV are strong enough to generate shape differences in the mouse isocortex and cerebellum, whereas, delipidation and PFA induces tissue shrinkage, especially in the Isocortex, cerebellum, and hindbrain regions (Fig. 2). Those differences make the building of a dissected, cleared tissue-specific atlas necessary.

Vascular development

The TubeMap pipeline developed by [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF], was adapted to the early developmental stages. From P1 to P9, brains were embedded in agarose cubes in order to fit the light sheet sample holder. The brains were transparized and then registered to the previously developed atlas. The resulting graphs were then characterized.

(Fig. 3) The amount of vessels increases in waves. The first one takes place between P3 and P7 and the second one between P14 and P21. As a redundant network, the evolution of the loop's shape is also a good feature to monitor. At P1 most of the loops are very small. During successive growth waves, the number of loops increases and their sizes increase with tissue expansion. The distribution of the loop's length becomes bimodal. At the end of the development, the whole distribution is shifted to a monomodal large loop distribution.

Vessels tip points

The Blaze light sheet microscope has enough signal resolution to show the filopodia tip cell at the end of the 1 degree node of the vessels, allowing us to label them as sprouting vessels. Measure based on the radius, length, and vessel intensity signal allows us to fully filter them out in order not to introduce topological aberration to our graphs.(Fig. 4). By modulating the clipping value in the TubeMap pipeline we can easily keep or remove the filopodia from the dataset.

Regional development dynamics

Voxelization of the graph branch points and end points shows that different regions have different time windows for development. By looking at the sprouting vessel density over time, we were able to classify the different regions according to their growth dynamics. Regions belonging to the same somatosensory network seem to co-develop synchronously.(Fig. 5 and Fig. 6). Therefore the synchronized development of these systems might correlate with the maturation of the sensory system it is related to.

Arterial development

The arterial network is also immature at birth and shows densification over the developmental phases.(Fig. 7) Until P14, the cortex shows a lack of penetrating arteries compared to the other meta regions. At P14 it looks like the arterioles and venules are overrepresented, and later on (1 month old) cleaned up only leaves the "useful" branches.

Orientation settings

In [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF] cortical vessels are shown to display preferential orientations depending on the brain region and cortical layer. The typical cortical organization is : an over-representation of penetrating vessels in layer 2/3, and an isotropic vessel orientation distribution in layer 4. In lower layers, an enrichment in planar vessels is observed. In the cortex, orientation is described relative to the surface since the region can be locally approximated with a 2D geometry. However, this concept could hardly be translated to the other brain metaregions where the surface displays more circumvolutions. It is interesting to observe that in the cortex, the main radial direction correlates with the penetrating artery's orientation. Therefore another, more generalizable way to compute vessels' orientation is to use the local arteries' and veins' direction as the normal vector from which the angle of the vessel will be compared.

To do so we compute the driving vector of every arterial or venous segment and use it to interpolate in any coordinates what would be the main flow orientation. Each vessel's orientation is defined by computing the angle it makes with the flow direction. (Fig. 8) During development, this layered cortical organization only starts to set place between P9 and P12. Before that, the network seemed to be mostly composed of penetrating vessels and immature arteries/veins couples. In most cortical regions, the complexification of the vascular net comes in two successive waves. The first one around P5 sees the sprouting of many "in-flow" vessels aligned along the arterial/venous flow direction, whereas the second wave, much more massive around P14 can be characterized by the isotropic sprouting of capillaries across the full cortex.

The reason behind the specificities of these orientations remains an open question. As the distribution of the flow and velocity values do not change much layer-wise, it can be hypothesized that it is needed to help redistribute the blood efficiently toward specific regions when needed.

The origins of this patterning are also unknown. Since it does not seem to be affected much by changes in neuronal activity (Fig. 14), the hypoxic theory is not likely to be valid. However, it could be due to some molecular gradient set in place during the development of the different brains structure. The rise of myelination, astrocytic co-development, pericytes maturation, thalamocortical wiring, sensory inputs, and regional differentiation are as many candidates potentially playing a role in the final organization of the vasculature.

Vessels communities

The idea of finding communities in a network is not new. In a natural network often some nodes are more similar to one another compared with its other components. In vasculature in particular it makes sense to consider that it can be subdivided into smaller arterio/venous domains. The idea of graph clustering or community partitioning is to divide large scale networks into sub networks, by maximizing a graph feature representing the ratio of intra versus inter communities edges. Here, we aim to understand the vascular organization. The previously described centrality measures can help , as graph features, to predict how fluid might be directed through this blood transportation network, and identify the key structure allowing precise distribution and regulation.

For instance, the vascular network is believed to be subdivided into veinal vs arterial domains, or simply subnetwork irrigated mostly from one of the penetrating arteries, where veins are just frontier-based «waste» collectors [START_REF] Kennel | From whole-organ imaging to in-silico blood flow modeling: A new multi-scale network analysis for revisiting tissue functional anatomy[END_REF] (Fig. 9) Because of the progressive complexification of the capillary bed, and the increased number of arteries and veins, the partitioning of the vasculature through development is susceptible to undergo deep changes with time.

Changes in the modularity and the number of clusters might be a good indicator of the resiliency and robustness of the network and could help to detect critical periods of the development process. In (Fig. 9) we show that the modularity of the network remains quite stable over time whereas the number of clusters increases. However, there are some abnormal timepoints, displaying discontinuities in the growth: in particular P9 (thalamus and striatum) and P21 (cerebellum and hindbrain) that correspond to inflection points on their network density during development.

DISCUSSION

We present a strategy to generate and analyze datasets of the developing mouse brain. With the multiplication of the clearing technique application, There has been a recent surge of interest in understanding the mechanism behind vascular and neuronal networks co development (Coelho Santos et al 2021). We show that the previously developed pipeline is flexible enough to be adapted to other structures than the only adult vasculature. We produce a highly granular and collaborative atlas of at 8 developing timepoints, and characterized the correctness of iDISCO clearing on tissues.

To address the diffusion anisotropy from [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF], we validated the use of conjugated antibodies to increase the signal on noise ratio in the center of the tissues (Sup Fig 1). Another issue to tackle concerns vessel filling. Since the antibodies CD31, Podocalyxin and SM22 on target proteins located on the surface of the vessel walls. Therefore, for large vessels, the signal fades in the centerline of the vessels. This results in empty vessels in the binarization, and may lead to topological aberration after the skeletonization part of the TubeMap pipeline.

Consequently, a DNN was trained in [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF] to detect those empty vessels in the binary mask and to fill them up. In the brain, the vast majority (90 %) of vessels has a radius between 2 to 5 pixels. Then most arteries and veins(roughly 10 % of the total vessels), forming the "large" brain vessels set, have a radius between 10 to 25 pixels. Finally some exceptionally surface veins and arteries (less than 1 % of the vessels) have a radius over 30 pixels. Because the DNN was trained on vascular graphs generated with blood still inside the vessels, the amount of extra large vessels in the training dataset was underrepresented even though they are the most susceptible to need to undergo filling. Therefore, even though the vast majority of empty vessels were processed in [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF], improvement could be made in order to fully catch the x-large surface vessels such as the basilar artery for instance.

MATERIALS AND METHODS l perfusion and tissue processing

All animals were sacrificed by Pentobarbital overdose (200mg/100 g of animal). Intracardiac perfusion was then performed with a peristaltic pump (Gilson, USA). Blood was washed by infusing 30mL of cold PBS followed by tissue fixation with 30mL of cold 4% Paraformaldehyde (Electron Microscopy Sciences, USA) diluted in PBS.

Then the brains were quickly dissected taking special care to preserve intact the structure. Finally, the brains were postfixed for 2 hours at room temperature by immersion in 4% para-formaldehyde, and stored in PBS at 4 C until further processing.

l atlas the autofluorescence of each sample is acquired in the light sheet microscope with a voxel resolution of 5 um × 5 um ×5.9 um. The autofluorescence is then downsampled to the allen atlas isotropic resolution (25 um × 25 um × 25 um) for the sake of coherence.

l alignement

We used an autofluorescence template from the ABA at 25mm generated from serial 2-photon tomography. To achieve the alignment, TubeMap first resamples the auto-fluorescent image of the sample to the resolution of the atlas reference and then aligns the image in 3d non-linearly. TubeMap integrates the Elastix package functionality [START_REF] Klein | elastix: a toolbox for intensity-based medical image registration[END_REF] for the non-linear alignment. For the vasculature data a hierarchical estimation of b-spline transformations between the two reference images is used together with a cross-entropy measure to quantify the local alignment quality. For the vasculature data, the reference image is acquired in tiles directly after each tile of the other channels and thus the same alignment and stitching layout can be used to assemble the full reference image. To correct for small misalignments between reference and data images a rigid transformation between both is estimated using also a cross-entropy measure. The resulting joint transformation from data to sample reference to atlas reference is used to transform the coordinates of the graph points onto the reference frame of the atlas. In the final step, the transformed positions of the graph are used to extract the atlas annotation and any other atlas information.

l template creation

To briefly give a quick overview of the process:

1 full brain autofluorescence acquisition 2 selection of the best brain (centered and symmetric) 3 only keep the right hemisphere of all samples 4 mirror the left hemisphere of all samples 5 aligned all the hemispheres on the selected one using the elastic transformation framework within Tubemap.

To create the atlas at each development timepoint, we acquired roughly 10 brains for each timepoint. Both hemispheres were acquired. Since mice brains are considered not to be lateralized, we used the right hemisphere as a template for annotation , and mirrored the opposite hemisphere in order to add samples to create the template by averaging them together. autofluorescence acquisition all autofluorescence of all timepoints were acquired the same way, in dorso-ventral orientation in the Blaze light sheet microscope from LaVision, at resolution 5×5×5.9 um/pixel, with both hemispheres.

alignment autofluo

We use the two steps alignment pipeline provided within Tubemap. First a rigid transformation is used to roughly overlap the autofluorescence on the template. Then a b-spline transformation between the two images is used for non-rigid deformations.

averaging strategies

After having align all the autofluorescence images together, they are average, the idea behind is that all small imperfections (tissue damages…) and brain inter individual dissimilarities are going to be blurred out by averaging them. All images are normalized in intensity using classical normalization, and intensity histogram matching. They are then stacked together and the resulting template is the median image. smoothing Since the coronal view was used to draw the regions, their surface presents a jerky aspect along the orthogonal axis (dorso-ventral, and medio-lateral). To smooth out the regions isotropically, each region was extracted based on its label. then a 3D mesh was created and smoothed using the laplacian smoothing algorithm provided by the PyVisa python library. the principle is that each vertex position becomes the average of its neighboring vertices, which has the effect of smoothing the mesh surface. Following the TubeMap framework, each brain to align should have an autofluorescence image downsampled to the atlas resolution , and the downsampled version of its stitched binary version of its vasculature. first the rigid transformation between both autofluorescence and downsampled binary files are computed,and then the elastic transformation between the autofluorescence and the template. The joint transformation is then applied to map the graph vertex coordinates into the atlas file and extract the regions information.

l statistics

Most statistics are computed with the Scipy library (https://scipy.org) and Numpy library. Some specific features from the graph-tool library (https://graph-tool.skewed.de) are also used to compute graph specific metrics. data are represented as mean and standard deviation. Mann-Whitney tests are used to assess if the mean difference between two groups is significant or not.

loops extraction

Given the redundant aspect of the vascular graph, loops are one crucial aspect to take into account when analyzing the graphs. but because of its redundancy, all parts of the vasculature is more or less part of a loop should it be large enough. Therefore we need to set an upper threshold in terms of size to decide what we consider a loop in the graph or not.

Based on total length criteria, we define as loops any circular subgraph composed of 3 to 7 vertices. A small collection of loops of 3 to seven vertices is created. For each number of vertices, we used the graph-tool isomorphic search function to obtain all the corresponding isomorphisms in the graph. The implementation is based on Cordella et al 2004Cordella et al (doi:10.1109Cordella et al /tpami.2004.75) .75) implementation that scales well to large graphs.

Branch point density, branch density, branch length and tortuosity

Branch point density and branch density were defined as the number of branch points or branches per volume in each brain region,respectively. The volume of each brain region was estimated by summing the volume of all voxels belonging to that region using the reference atlas annotation. The length of each branch was measured as the sum of the Euclidean distances between consecutive pairs of coordinates that define the branch as extracted from the raw graph. Branch tortuosity was defined as the ratio between the length of the branch divided by the Euclidean distance between its endpoints. total length of arteries because the arterial graph is disconnected and noisy, branching point density is not well adapted. A less prone to error measure is the total length.

density

The volume used to compute the density is computed by counting the number of voxels corresponding to any regions according to the atlas it is aligned with. then this number is multiplied by its resolution in order to get the volume expressed in um^3. voxelization It refers to the heatmap representation of the spatial graph. For instance, for a branch point voxelization a sphere of radius r is drawn around the atlas coordinates of each branch point, with a value 1. The spheres add to each other. The local value of the pixel is the number of spheres that overlap on top of this voxel, thus scoring locally the density of branch point.

l orientation

In Kirst et al 2020 we saw that vessels show specific organization, especially on the cortex. These specificities make the orientation of the vessels a key feature to compute.

local normal

The first method aims to compute the orientation on the cortex. Because of its relatively planar aspect the cortex can be locally approximated as a flat surface. We extract each cortical subregion, and approximate the upper layer of it as a surface by fitting a plane equation to the point cloud generated by its vertices. the normal to this surface is computed.

Alternative methods to increase the accuracy of the normal vector are:

To compute the normal of its lower layer and take the average of the two vectors as the effective normal.

To use the vertex information property containing the distance to the surface to compute the direction of distance gradient and adjust with the orthogonal to the surface.

To compute the main arteries orientation which are known to be penetrating in the cortex.

arteries based

The limitation of the previous method lies in the fact that it only works for the cortex, but cannot be generalized to the rest of the brain. Since the arterial direction matches with the normal to the surface on the cortex, we can use them as a reference to estimate the orientation of each capillaries. They can be either collinear to the arterial flow, or orthogonal to it. the methods is the following: First we compute the orientation vector for each artery labeled edges, which define the blood flow vector field.

Then a grid is created that covers the surface of the brain. Based on the previously computed partial vector field, The 1D value in X, Y and Z is estimated using interpolation function F. Using F, the hypothetical blood flow direction is computed on each capillary labeled edges. it corresponds to the normal vector from with the orientation is to be computed. For each capillary we compute the angle between the vessel and its normal.

Every vessel with an angle inferior to 40 degrees is considered collinear to the blood flow or "in-flow", and orthogonal or "cross-flow" otherwise.

flow interpolation

Similarly to the previously described method:

First we compute the orientation vector for each artery labeled edges, which define the blood flow vector field.

Then a grid is created that covers the surface of the brain. Based on the previously computed partial vector field, The 1D value in X, Y and Z is estimated using the interpolation function F. Using F, the hypothetical blood flow direction is computed on each point of the grid. l clustering hierarchical region dynamics clustering Clusters are computed using the ward clustering algorithm (Ward et al. 1963). The method aims to minimize the inertia of the system by iteratively clustering together the pair of elements or clusters that minimizes the increase of the total inertia.

57

V. Discussion

In this dissertation, I present a methodology to generate and analyze graphs of the adult and developing mouse brain vasculature. I applied this pipeline to describe the organization of the adult cerebral vasculature, and how its remodeled following a stroke, a viral infection, or a loss of sensory inputs. While the current pipeline offers countless opportunities to better understand the organization of the normal and pathological cerebral vascular network, there is still a lot of development to do to improve our capabilities to run such analyses.

On one hand, the quality of the segmentation itself could be improved for the very large vessels. Moreover, the segmentation of the surface pial vasculature also represents a significant hurdle, not addressed in the current pipeline.

On the other hand, most of the analysis done on vascular graphs are still simple, and do not yet harness all the information present in these reconstructions. Novel methodologies to analyze the brain vasculature, adapted to the very large size of the data, should be developed to extract more features of its organization, and determine which aspects of this organization are changed in pathological conditions. To improve these analyses, the quality of the data produced must also be improved, for the molecular annotation of veins and arteries.

In the first part of this discussion, I will present several ideas to improve technically the ClearMap segmentation.

In the second part, I will present new mathematical tools developed to study vascular plasticity. Finally, I will discuss how preliminary data obtained at the end of my PhD could illuminate the life-long interactions between neuronal activity and blood vessels.

V.I Improve the quality of vascular reconstructions.

The possible areas of improvement for our pipeline could be divided in 4 major directions:

1. Improving the segmentation of small vessels in the center of the brain 2. Improving the reconstruction of very large empty vessels 3. Improving the reconstruction of surface vessels

Improving the venous and arterial discrimination

Improving the segmentation of small vessels in the center of the brain: Improving the occasional failure to segment deep and small capillaries would require solving the diffusion gradient we observe from the immunostaining. Indeed, the signal originating from deep vessel is weak and can lead to interruptions in the reconstruction. To address the diffusion anisotropy from [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF], we validated the use of conjugated antibodies to increase the signal-on-noise ratio in the center of the tissues (Fig 16). Removing the secondary staining step made the diffusion gradient flatter in brains processed with directly conjugated antibodies. While this resulted in a slightly dimmer signal at the surface of the brain compared to the traditional dual-stage labeling approach, the signal in the center of the brain was brighter than when using a primary-secondary conjugated pair. This enabled a more reliable segmentation of the vasculature in the thalamus than in our initial datasets. Other strategies could be employed, for instance the promising use of tail-vein injected fluorescent lectins, that circulate in the blood stream and label the inner wall of all vessels. This strategy would in principle lead to a perfect labeling throughout the depth of the tissue but would have a potential issue for the final reconstruction. Indeed, capillaries can be constricted and closed for either physiological or pathological reasons, and the circulating dye may not diffuse effectively in the whole vascular arbor. Therefore, injected dye-based methods, while providing a very signal-to-noise ratio, may not reveal the complete brain vasculature.

Improving the reconstruction of very large empty vessels: Another issue left to tackle on the segmentation parts of the pipeline concerns the reconstruction of very large vessels. Most of the large vessels (larger than about 20µm in our data) are located at the surface of the mouse brain, but some are present in the parenchyma.

The antibodies used for our vascular labeling (CD31, Podocalyxin, and SM22) target proteins that are located on the inner surface of the large vessel walls (the endothelial cells). Therefore, for large vessels, the center is empty of signal. These vessels therefore appear as tubes in the dataset. This results in empty vessels in the binarization and may lead to topological aberration after the skeletonization part of the TubeMap pipeline. In our published work, I trained a DNN [START_REF] Kirst | Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature[END_REF] to detect those empty vessels in the binary mask and to fill them. In the brain, most vessels have a radius between 2 to 5 pixels. Most arteries and veins (about 10 % of the total vessels), forming the "large" brain vessels set, have a radius between 10 to 25 pixels. Finally, very large veins (less than 1 % of the vessels) have a radius over 30 pixels. I trained the DNN on vascular graphs generated with blood still inside the vessels to improve their segmentation (hence corresponding to a dataset more like those generated by dye injections). The number of large vessels in the training dataset is underrepresented owing to their natural sparseness relative to capillaries. Most medium-sized empty vessels are correctly filled by the initial network I trained (corresponding to venules and arteries around a size of 10 to 20µm in diameter). However, as the size of the vessels increases, their detection and filling by the DNN decrease. This leads to the presence of poor reconstructions of large vessels in our data. Since the wall of these vessels is still segmented, but when they are not recognized by the DNN, their center stays hollow, the topological thinning algorithm fail to reconstruct a center line and loops around the hollow center, creating structures akin to a fish net in the final graph. Improvements should therefore be made to fully catch the larger veins and surface arteries, such as the basilar artery for instance. An alternative solution to the "vessel filling" network, would be to detect the poorly reconstructed vessels directly from the graph. The topological thinning of a hollow vessel indeed generates a very peculiar structure (the "fish net") that is easy to identify visually when plotting in 3D the vascular graph. It should be therefore possible to design a neural network to detect these occurrences and correct the graph structure to reconstruct a continuous vessel. Such neural network may run faster on the graph data than a convolutional U-net running on the very large image data.

Improving the reconstruction of surface vessels: During my PhD, I focused on the study of blood vessels located in the brain parenchyma, excluding the circulation network present at the surface of the brain. This is partly because of the difficulty to reconstruct the surface pial network with our current tools. Unlike the parenchymal circulation, which is organized in 3D loop structures, the surface pial circulation is organized on the 2D brain surface. This implicates that surface vessels are often running on top of each other's: this could be a small arteriole running along a larger vein, or even a large artery crossing over a large vein. There are many factors that limited my capacity to correctly segment the vasculature on the brain surface. The first is that the resolution of the light sheet microscope used, 1.6µm/pixel lateral, and 2.5µm/pixel axial, is insufficient to discriminate the crossings of surface vessels, especially when vessels are running along each other's.

Upgrading the imaging from 4X to 12X (hence 0.5µm/pixel lateral) may improve this issue, but this wouldn't change the axial resolution. Scanning a brain at 12X used to be a daunting endeavor due to the very long scan-times it would entail. However, a recent development on how the motor sequenced is programmed into the light sheet microscope has been upgraded, leading to a over 10-fold decrease in scan times. This may enable the use of higher magnification objectives soon, necessary to improve the surface reconstructions.

Another point to consider is that the strategy used for the centerline reconstruction relies on a 3D topological thinning, which would produce artefactual results in case of vessels crossing or running close. Indeed, it would only take a touching pixel between 2 vessels to generate an aberrant crossing with this strategy. I started In Figure 18, I measured the modularity and number of modules in the main brain regions across post-natal development, from P1 to P30. The high variance observed in the modularity of certain regions at certain time points may be due to the preliminary nature of these data and may be resolved by adding higher number of cases with better reconstructions. Nevertheless, we can appreciate here that the modularity of the vascular network stays stable during development across all regions. The number of modules increases over time as expected. This suggests, if this trend is confirmed with more data points, that the modularity may not be a measure adapted to understand variations during normal brain development, as the system seems to enforce a stable modularity in normal conditions.

I also tested whether the measure of graph modularity could be relevant in a situation where we already had measured a natural vessel loss, in the auditory cortex of the deaf Otoferlin mutants. The vessel loss in the auditory system of the Otoferlin mutants always fluctuated between 10% and 20% of loss in the mutant compared to littermate controls. The question was whether the lower number of vessels in that case would lead to a change in the organization of the vascular modules in those regions. I had 2 datasets I could directly use, one captured at 1 month, and one at 6 months.

A major unanswered question from the initial analysis I had run based on simple metrics is whether the vessel loss affects different clusters of the whole network differently, or whether the loss is equally spread over the whole region. To test this, I am proposing to use a simple rerouting model that randomly removes edges while preserving the biological characteristics of the original vascular network (Fig. 19). We applied this strategy to the vascular network of the auditory cortex, for which I have experimental data points in normal and deaf animals. The idea of this strategy was to check if a random edge removal could capture the differences in vascular organization of the deaf mice. To run this comparison, I used the SBM methodology to measure the modularity and number of modules in both the artificially randomly pruned graph (from wild type littermates), and the mutant graphs.

The random pruning strategy was able to reproduce the network properties of the mutant mice. After removing around 12% of the edges from the healthy mice (Fig. 19), I obtained on average similar numbers of links, modules, and modularity values as observed in the mutant mice (Fig. 19).

In sum, this work represents a preliminary step toward the identification of useful measures to describe the network mechanisms of vascular plasticity. In the initial tests presented here, it looks like the use of SBM to estimate the modularity of the vascular network (which is different from the actual vascular domains, as defined by the arterio-venous flow), is only useful in certain cases where a pathological condition affects the normal structure of the network. Interestingly, this may indicate that manipulating neuronal activity may impact the modularity and number of modules of a brain region. If this could be reproduced in different models, this could be an interesting method to evaluate vascular plasticity, as it is relatively light computationally, and very robust to potential artefacts present in the data, unlike flow measures. Moreover, it doesn't require to have labeled veins and arteries, which, as I pointed out above, is still challenging. vesicle transportation in the auditory sensory hair cells, resulting in total hearing loss. This induces a lack of sensory inputs in the auditory system during the mouse development and allowed me to measure its effect on the vascular network in the whole brain, both in the auditory system and elsewhere.

V.3.1 Effect of neuronal activity on vascular development

My repeated analysis of these mutants showed a significant decrease of roughly 12 % of the number of branch points in the auditory cortex, in the lower layers, on the one hand, but also in subcortical auditory relays such as Inferior Colliculus, lateral lemniscus of the Thalamus, and the Olivary complex (Fig 21).

Surprisingly, I found that when considering the orientation of blood vessels, the proportion of in-flow and crossflow vessels in the different cortical layers remained unaffected, suggesting that neuronal activity is involved in the densification process of the vasculature, but not its intrinsic organization structure. It could suggest that the structure of the vascular network, or at least the orientation of the vessels, is regionally patterned by molecular cues or gradients, while the number of vessels could be controlled in a more adaptive way.

Another interesting insight from these data was to notice that, on top of the decrease in the vascular densities in the auditory regions, there was a transient increase in the somatosensory regions during early development, which disappeared later. A lack of auditory inputs is known to potentiate the activity of the somatosensory cortex, as both are strongly connected. At P14, I couldn't measure a significant decrease in the vascularization of the auditory system, but the somatosensory system was already hyper-vascularized. Because this

V.3.3 Effect of late deprivation and aging on the vascular system

As the developing brain is known to be plastic, it remains unclear whether the same applies to the adult brain.

Unfortunately, the experiments testing the effects of late deprivation, both in the auditory system and in the somatosensory system didn't yield interpretable results, as the variability was too high.

However, it was interesting to note that, in the early-deprived otoferlin mice, the degradation of the auditory system vasculature, but also other systems, was still an ongoing phenomenon as the mice aged, with vascular reconstructions going all the way through 10 months-old animals. Indeed, we measured a worsening of the vascular density in many brain regions connected to the auditory system, such as the Ventral Tegmental Area, the pre-frontal cortex, or associative cortices. Strikingly, the only region spared from this degradation seemed to be the somatosensory system, which may be more active during the adult life of the deaf mice, as our Fos data suggested (Figure 23).

It would be important to repeat experiments of late deprivation, in models where neuronal activity is continuously depressed, to establish measures of vascular density across time and aging. Finding the right models of adult activity deprivation may be challenging: optogenetic or chemogenetic manipulations are not permanent and require surgical interventions. Kir2.1 expression in neurons to depress activity permanently could be an avenue to explore. Moreover, the use of transcranial noninvasive optogenetic implants may be a way to circumvent the pitfalls of traditional optogenetic stimulations via optic fibers.

While the exact effect of aging on the vascular structure is yet to be uncovered, the preliminary data presented here suggest the presence of life-long adjustments to the blood vessel density in the brain. This could reshape the way we think about how aging or brain disorders such as depression may affect the metabolic supply to neurons.

network. While the network retains a high level of plasticity during development and in pathological conditions, a crucial piece of the puzzle remains: how adult brain adaptability and aging affect this vascular network.

Understanding this relationship is vital, and cutting-edge technologies that enable 3D reconstructions of the vascular network are at the forefront of this effort. Adult plasticity, the brain's capacity to reorganize and adapt, may have close ties to vascular adaptations and changes. This could support processes like learning, memory, and recovery from injuries by tuning the metabolic support to neurons in a finer way than with pure vasomotor control.

Aging poses its own set of challenges. The brain undergoes structural and functional changes as it ages, some of which may be directly linked to the blood vessels. However, reports of vascular changes as the brain ages are still not always in agreement over whether only pericytes, or other vascular cells, are affected. These agerelated shifts, including declining vascular health and reduced adaptability, raise concerns about cognitive decline and be a major factor driving an increased risk of neurodegenerative diseases.

This is where 3D vascular network reconstructions come into play. They offer a powerful tool to visualize and study the intricate changes occurring within the cerebral vascular network as the brain adapts and ages. These technologies provide us with a clearer picture of how the network evolves over time, but still require much more work in terms of accuracy, and the downstream mathematical tools to take advantage of them.

In conclusion, the cerebral vascular network is far from static; it's a dynamic part of the brain that evolves with age and adapts to new challenges. To grasp the full picture of how adult plasticity and aging impact this network, we must still rely on cutting-edge 3D reconstructions, but improve the quality of these reconstructions.

This is because we need to quantify minute changes to the vascular network, which may have an outsized impact on brain health and function.

ANNEX

A versioning tool for collaborative atlas enrichment

As stated before in the second paper presented, one of the perks of building developmental atlases at so many timepoints is the novelty and the importance of the resources for future works on pups it will allow to unlock.

However due to the 3D nature of the datasets, the number of regions, and the difficulty of precisely manually segmenting them, it represents a long and tedious task. Moreover depending on the system you might want to study, the scale, the granularity, and the regions you might want to look at can be very diverse. All the regions do not exist in pups as they are in adults. And segment all the Allen brain regions on the pups' brain seems like an overkill and mostly irrelevant approach. Therefore we released an online collaborative tool based on versioning concepts just like github. Allowing each user to annotate regions of interest that were not already labeled and share their modification with the rest of the community.

What is versioning?

Versioning

Introduction

As of September 2020, SARS-CoV-2 has infected >25 million people globally. While a majority of COVID-19 patients present with respiratory symptoms, neurological involvement, including impaired consciousness and headache, have been reported in patients [START_REF] Mao | Neurologic Manifestations of Hospitalized Patients With Coronavirus Disease 2019 in Wuhan, China[END_REF]. To date, human autopsy studies have identified viral RNA transcripts in brain tissues [START_REF] Puelles | Multiorgan and Renal Tropism of SARS-CoV-2[END_REF][START_REF] Solomon | Neuropathological Features of Covid-19[END_REF] and viral proteins in the endothelial cells within the olfactory bulb [START_REF] Cantuti-Castelvetri | Neuropilin-1 facilitates SARS-CoV-2 cell entry and infectivity[END_REF] in people who succumbed to COVID-19. Several reports have come out indicating SARS-CoV-2 infection of the central nervous system (CNS) cells [START_REF] Bullen | Infectability of human BrainSphere neurons suggests neurotropism of SARS-CoV-2[END_REF][START_REF] Jacob | Human Pluripotent Stem Cell-Derived Neural Cells and Brain Organoids Reveal SARS-CoV-2 Neurotropism Predominates in Choroid Plexus Epithelium[END_REF][START_REF] Pellegrini | Long-Term Neurological Threats of COVID-19: A Call to Update the Thinking About the Outcomes of the Coronavirus Pandemic[END_REF][START_REF] Ramani | SARS-CoV-2 targets neurons of 3D human brain organoids[END_REF][START_REF] Yang | A Human Pluripotent Stem Cellbased Platform to Study SARS-CoV-2 Tropism and Model Virus Infection in Human Cells and Organoids[END_REF][START_REF] Zhang | Long-term bone and lung consequences associated with hospital-acquired severe acute respiratory syndrome: a 15-year follow-up from a prospective cohort study[END_REF]. However, there are still many unknowns regarding the frequency or the consequences of neuroinvasion.

Understanding the full extent of viral invasion is crucial to treating patients as we begin to try to figure out the long-term consequences of COVID-19, many of which are predicted to have possible CNS involvement [START_REF] De Felice | Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) and the Central Nervous System[END_REF][START_REF] Heneka | Immediate and long-term consequences of COVID-19 infections for the development of neurological disease[END_REF]Pereira, 2020;[START_REF] Zhang | Long-term bone and lung consequences associated with hospital-acquired severe acute respiratory syndrome: a 15-year follow-up from a prospective cohort study[END_REF]. Because the CNS is not the primary organ affected by SARS-CoV-2, studying neurological disease in COVID-19 patients systemically provides several challenges, including having only a subset of the population of patients with neuroinvasion, lacking technology to sample CNS tissues directly, and distinguishing direct neuroinvasion versus systemic viremia within the brain. Therefore, robust, reliable model systems are required to answer the questions underlying SARS-CoV-2 neuropathology. During the Zika virus (ZIKV) pandemic, several groups used human brain organoids to answer key questions regarding ZIKV neuroinvasion and its consequences [START_REF] Garcez | Zika virus impairs growth in human neurospheres and brain organoids[END_REF][START_REF] Qian | Brain-Region-Specific Organoids Using Mini-bioreactors for Modeling ZIKV Exposure[END_REF][START_REF] Xu | Zika virus infection induces RNAi-mediated antiviral immunity in human neural progenitors and brain organoids[END_REF]. Using a similarly well-characterized human brain organoid model [START_REF] Lancaster | Generation of cerebral organoids from human pluripotent stem cells[END_REF][START_REF] Lancaster | Cerebral organoids model human brain development and microcephaly[END_REF], we test the infection capacity of SARS-CoV-2 in CNS tissue. Using single-cell RNA sequencing (RNA-seq), we uncover the transcriptional changes caused by SARS-CoV-2 infection of neurons.

Beyond the neuroinvasive potential of SARS-CoV-2, the question remains whether ACE2 is the main route of entry of SARS-CoV-2 into neuronal cells and what strategies might block viral infection. ACE2 expression in the CNS, and in neurons in particular, is still unclear [START_REF] Xia | Angiotensin-converting enzyme 2 in the brain: properties and future directions[END_REF]. Moreover, in addition to ACE2, other cofactors such as TMPRSS2 [START_REF] Hoffmann | SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor[END_REF] and Neuropilin-1 (Cantuti-Castelvetri et al., 2020;[START_REF] Daly | Neuropilin-1 is a host factor for SARS-CoV-2 infection[END_REF] seem to affect infection rates, but it is unclear if these are also required for neuronal infection. Thus, we used transcriptional profiling along with blocking antibody studies to demonstrate the requirement of ACE2 and SARS-CoV-2 spike protein to infect neurons.

To gain in vivo relevance for these findings, we examine the neuroinvasive potential of SARS-CoV-2 using mouse models of SARS-CoV-2, and observe vascular remodeling in infected regions, independent of vascular infection, providing a valuable tool that can be used to dissect out the consequences of CNS infection. Finally, by examining postmortem COVID-19 patient brain tissues, we provide evidence of neuroinvasion by SARS-CoV-2 and identify associations between infection and ischemic infarcts in localized brain regions.

Results

Modeling SARS-CoV-2 neuroinvasion and cellular death using human brain organoids To dissect the mode and consequences of infection, we first established the neuroinvasive potential of SARS-CoV-2 in a human brain model system. We used human induced pluripotent stem cell (hiPSC) lines (Y1 and Y6), derived from healthy individuals (Fig. 1), to generate forebrain-specific human neural progenitor cells (hNPCs). In culture, we observed replication of SARS-CoV-2 in 2-wk-old hNPCs, with peak viral titers as early as 12 h postinfection (hpi; Fig. S1, A andB). In addition, Tdt-mediated dUTP-biotin nick end labeling (TUNEL) staining indicated increased cell death (Fig. S1 A). Next, we generated hiPSC-derived brain organoids [START_REF] Amin | Building Models of Brain Disorders with Three-Dimensional Organoids[END_REF][START_REF] Pellegrini | Long-Term Neurological Threats of COVID-19: A Call to Update the Thinking About the Outcomes of the Coronavirus Pandemic[END_REF][START_REF] Qian | Brain-Region-Specific Organoids Using Mini-bioreactors for Modeling ZIKV Exposure[END_REF][START_REF] Velasco | Individual brain organoids reproducibly form cell diversity of the human cerebral cortex[END_REF] to model the SARS-CoV-2 infection of neuronal cells in 3D. We confirmed a dorsal cortical identity of organoids by immunostaining for markers FOXG1, PAX6, and CTIP2 (Fig. S2 C). Similar to a recent report [START_REF] Ramani | SARS-CoV-2 targets neurons of 3D human brain organoids[END_REF], we observed infection of neuronal cells in 9-wk-old organoids as early as 24 hpi, with significantly increased number of SARS-CoV-2-positive cells at 96 hpi (Fig. 2, A-C;a n dFig. S1 C). Although the majority of the SARS-CoV-2-infected cells were localized within MAP2-positive cellular fields of mature neurons (Fig. 2 B [2] and Fig. S1 C), we also observed infection of SOX2-positive neural stem cells with bipolar morphology and cells localized around the neural tube-like structures (Fig. 2 B [1]). We observed an increase in SARS-CoV-2-positive cells 96 hpi compared with 24 hpi (Fig. 2 C and Fig. S1, C-E). 96 hpi, we observed a widespread infection, mostly limited to the regions with high cortical cell density in the organoid (Fig. 2 D). Using electron microscopy, we visualized viral particles within the organoid (Fig. 2 E), with discrete regions of high-density virus accumulation (Fig. 2 E [1]) and other regions showing viral budding from ER-like structures (Fig. 2 E [2][3][4]), suggesting that SARS-CoV-2 is able to use the neuronal cell machinery to replicate. Organoid infection resulted in extensive neuronal cell death; strikingly, however, the majority of TUNEL-positive cells were SARS-CoV-2 negative (Fig. 2, F and G; and Fig. S1 F), and only ∼15% of the cells infected with SARS-CoV-2 were TUNEL positive (Fig. 2 H). Increased cell death was correlated with a higher density of SARS-CoV-2-positive cells (Fig. 2 I), which, however, were not overtly dying (Fig. 2, G andH). Even within a single plane, we noticed the high-density SARS-CoV-2 area to have more TUNEL-positive cells (Fig. 2 J, yellow box) compared with low-density SARS-CoV-2 regions (Fig. 2 J, white box). Together, these data indicated that SARS-CoV-2 can infect cells of neural origin and suggested that infected cells can promote death of nearby cells. population changes during infection (Fig. 4). With an added SARS-CoV-2 annotation, SARS-CoV-2 transcript reads were localized to a variety of cell clusters (Fig. 3 A), demonstrating the widespread infectivity of SARS-CoV-2 in neurons, radial glia, and neuronal progenitor cells. Several cell clusters showed large changes in their representation within the organoid, such as cluster 1 (from 5% to ∼20% of the population), or cluster 7 (from ∼10% to nearly 0%; Fig. 4). However, this was not uniformly seen with all infected clusters, as cells in cluster 11 showed high infection rates (Fig. 3 G)b u tl i t t l ec h a n g ei n population representation (Fig. 4). This is consistent with findings from TUNEL staining (Fig. 2, E-G), which demonstrated minimal overlap between SARS-CoV-2-infected cells and those undergoing cell death. Alternatively, it is possible that infection resulted in differentiation or change in cell's phenotype in populations such as cluster 7, which is characterized to be a transitional cell type.

Next, we compared the changes in the transcriptome of SARS-CoV-2-infected organoids with another well-studied neurotropic virus: ZIKV. Differentially expressed genes (DEGs) after a brain organoid infection with ZIKV [START_REF] Watanabe | Self-Organized Cerebral Organoids with Human-Specific Features Predict Effective Drugs to Combat Zika Virus Infection[END_REF] showed almost no overlap with DEGs from SARS-CoV-2infected organoids (Fig. 5 B). In both cases, the transcriptome showed evidence of viral infection and invasion (Fig. 5 D). However, unique processes were enriched in each of the infections. SARS-CoV-2-infected brain organoid up-regulated pathways related to cell division, organelle fission, and metabolic processes, while ZIKV showed enrichment in type I IFN pathways (Fig. 5 D). SARS-CoV-2 induced a unique transcriptional state within neurons compared with ZIKV, consistent with the fact that SARS-CoV-2 induces a moderate IFNstimulated gene response in other tissues [START_REF] Blanco-Melo | Imbalanced Host Response to SARS-CoV-2 Drives Development of COVID-19[END_REF] and previous reports of specific virus replication being controlled by alternative pathways by neurons [START_REF] Daniels | The Nucleotide Sensor ZBP1 and Kinase RIPK3 Induce the Enzyme IRG1 to Promote an Song et al[END_REF][START_REF] Yordy | A neuronspecific role for autophagy in antiviral defense against herpes simplex virus[END_REF].

To understand the impact of SARS-CoV-2 infection on the brain organoid at the cellular level, we characterized cells that were infected with SARS-CoV-2 versus neighboring uninfected cells without the presence of SARS-CoV-2 transcript, in cluster 11, which comprised the highest number of infected cells 96 hpi. SARS-CoV-2-positive cells showed enrichment of genes corresponding to viral transcription, along with enrichment for metabolic processes including electron transport-coupled proton transport, cytochrome c to oxygen, and NADH to ubiquinone (Fig. 5 E Host ACE2 receptor is required for infection of neurons One of the ongoing questions regarding SARS-CoV-2 neuroinvasion is that mRNA levels of ACE2 appear to be very low in the CNS [START_REF] Li | Expression of the SARS-CoV-2 cell receptor gene ACE2 in a wide variety of human tissues[END_REF][START_REF] Qi | Single cell RNA sequencing of 13 human tissues identify cell types and receptors of human coronaviruses[END_REF][START_REF] Sungnak | SARS-CoV-2 entry factors are highly expressed in nasal epithelial cells together with innate immune genes[END_REF]. Indeed, our single-cell RNA-seq dataset demonstrated low levels of ACE2; it was, however, detectable in many clusters (Fig. S2 A).

In addition, we did not observe a correlation between the percentage of cells infected in each cluster to ACE2, TMPRSS2, or Neuropilin-1 expression (Fig. S2 B). However, it remained possible that ACE2 protein may be expressed on the cell surface to promote viral entry. Consistent with this idea, we found widespread expression of ACE2 protein in both MAP2-positive neurons and cells in the neural tube-like structures of the organoids (Fig. 6 A), indicating that the mRNA level of ACE2 does not accurately reflect ACE2 protein expression. In addition, we used postmortem human brain tissue to stain for neurons and ACE2 and found that neurons in the cortical gray matter colocalized with ACE2 staining (Fig. S2 C, yellow arrow), and we found other cells in the vicinity that also stained for ACE2 (Fig. S2 C, white arrow). To test the requirement of ACE2 for SARS-CoV-2 infection, we incubated organoids with an anti-ACE2 blocking monoclonal antibody before infection with SARS- CoV-2. We detected significant inhibition of SARS-CoV-2 infection upon pretreatment with ACE2 antibody compared with isotype control, indicating the requirement of ACE2 for infection of brain organoids (Fig. 6, B and F;an dFig. S2 C). Next, we investigated whether there are humoral antibody responses against SARS-CoV-2 in the CNS of infected patients, and whether antibodies present in the CNS can prevent infection of neurons. We analyzed cerebrospinal fluid (CSF) from a patient hospitalized with COVID-19 and acute encephalopathy (Fig. 6 C) and from a healthy control volunteer by performing ELISA against the spike protein of SARS-CoV-2 (Fig. 6 D). We detected IgG antibodies specific to the spike protein in the patient's CSF even at 100× dilution (Fig. 6 D). Using this patient's CSF, we performed a neutralization assay against SARS-CoV-2 and validated its efficacy in preventing brain organoid infection (Fig. 6 E). CSF-containing antiviral antibodies blocked SARS-CoV-2 infection in organoids (Fig. 6, E and F;a n dFig. S2 D).

Mouse models of COVID-19 confirm the neuroinvasive potential of SARS-CoV-2

To examine the consequences of SARS-CoV-2 infection in a more physiologically complete system, we examined the neuroinvasive potential of SARS-CoV-2 in vivo by using transgenic mice expressing human ACE2 under the K18 promoter (K18-hACE2; [START_REF] Mccray | Lethal infection of K18-hACE2 mice infected with severe acute respiratory syndrome coronavirus[END_REF]. Similar to previous reports of SARS-CoV showing neurotropism [START_REF] Mccray | Lethal infection of K18-hACE2 mice infected with severe acute respiratory syndrome coronavirus[END_REF][START_REF] Netland | Severe acute respiratory syndrome coronavirus infection causes neuronal death in the absence of encephalitis in mice transgenic for human ACE2[END_REF], we observed increasing viral titers in the brain of mice after intranasal administration of SARS-CoV-2 (Fig. 7, A andB). We next analyzed the distribution of the virus in the whole brain by immunolabeling against the nucleocapsid protein, clearing, and light sheet microscopy imaging using iDISCO+ (Renier et al., 2014;Fig. 7 and Video 1). 7 d after infection, the virus was widely present in neural cells throughout the forebrain (Fig. 7 C). The cortex was unevenly infected, as the infected cells were visible in columnar patches and in sensory regions, while the layer 4 was mostly devoid of infection (Fig. 7, Ca n dD ). We mapped the density of infected cells using Clear-Map (Renier et al., 2016) and confirmed that most brain regions contained a high density of infected cells, with the notable exception of the cerebellum (Fig. 7 E). Other regions also contained a relatively low density of infected cells, for instance the dentate gyrus, the globus pallidus, and cortical layer 4. Of note, wholebrain analysis of the virus distribution did not detect the presence of the virus in the vascular endothelium. To explore the possibility that the viral expression by neural cells could indirectly affect the organization of the vascular network, we double-labeled brains for the nucleocapsid protein and the vascular endothelium with CD31 and Podocalyxin, and used ClearMap again to reconstruct the whole-brain vascular network 7 d after infection (Fig. 7 F; [START_REF] Kirst | Mapping whole brain structure and activity at cellular resolution: iDISCO+ and ClearMap[END_REF]. We focused on the organization of the cortical vasculature. In infected brains, we focused on hot spots of infections represented as the number of detected infected cells and measured the density and orientation of all blood vessels. In many instances, the expression of the virus coincided with a disruption in the normal vascular topology expected in the cortex, with an important loss of the normal enrichment in radially oriented vessels in upper layers.

To dissect the consequences of SARS-CoV-2 infection in the CNS versus the respiratory system, we directed expression of human ACE2 using an adeno-associated virus vector (AAV-hACE2; [START_REF] Israelow | Mouse model of SARS-CoV-2 reveals inflammatory role of type I interferon signaling[END_REF] to the lungs via intratracheal delivery or to the brain via intracisternal delivery (Fig. 7 H). Mice were then infected with SARS-CoV-2 either intranasally or intraventricularly. Intranasally infected mice showed signs of lung pathology [START_REF] Israelow | Mouse model of SARS-CoV-2 reveals inflammatory role of type I interferon signaling[END_REF] but no weight loss or death (Fig. 7, I and J). However, intraventricular administration of SARS-CoV-2 resulted in weight loss and death, even at the challenge virus dose 100-fold lower than that used for intranasal infection (Fig. 7, I and J). Altogether, these results highlight the neuroreplicative potential and lethal consequences of SARS-CoV-2 CNS infection in mice.

Evidence of SARS-CoV-2 neuroinvasion in COVID-19 patient brain autopsies Finally, to determine whether there is evidence of SARS-CoV-2 CNS infection in COVID-19 patients, we examined formalinfixed, paraffin-embedded (FFPE) sections of three patients who died after suffering from severe COVID-19-related complications. All patients had been admitted to the intensive care unit, they were sedated and ventilated due to respiratory failure (days 3, 10, and 18 for patients 3, 1, and 2, respectively), and their difficulty to be weaned from mechanical ventilation indicated the severity and highly pathogenic nature of the disease course (Table 1). We sampled several regions of the brain and stained with H&E and for SARS-CoV-2 spike protein with a validated protocol and antibody we have used previously [START_REF] Hosier | SARS-CoV-2 infection of the placenta[END_REF]. While in control slides, there was no positive immunohistochemical staining for the SARS-CoV-2 spike protein, samples from COVID-19 patients showed distinct and specific spike protein staining, albeit to different degrees (Fig. 8 and Fig. S3). Neuropathological evaluations identified staining with antispike antibodies in cortical neurons (Fig. 8 A In patient 1, medium spiny neurons in the basal ganglia did not show staining, however, similar to substantia nigra neurons in the midbrain, which were also negative for the SARS-CoV-2 spike protein (Fig. S4). The cellular staining pattern showed diffuse cytoplasmic and perinuclear staining, along with concentrated regions within the cells, consistent with our findings using electron microscopy in the organoids, which demonstrated localized hotspots within the infected cell. Tissue samples from patients with cells that were lipofuscin laden or contained microinfarcts but without COVID-19 showed no reactivity to the SARS-CoV-2 antibody, signifying the specificity of the antibody used. Patient 3, who failed to regain consciousness after cessation of sedation, was diagnosed with severe global encephalopathy and was noted to have multiple diffuse microbleeds, as demonstrated by diagnostic imaging (magnetic resonance imaging). Upon histological examination, we found multiple microscopic ischemic infarcts in the subcortical white matter, ranging from acute to subacute, and with focal hemorrhagic conversion (Fig. 8 C and Fig. S5 A). These varying stages of infarction indicated a temporal sequence of continued ischemic events. Most infarcts showed signs of tissue damage and localized cell death, and positive viral staining was present predominantly around the edges of the infarct and to a lesser degree within the center. At the hyperacute stage of infarction, viral proteins were present in endothelium. These infarcts are present predominantly in subcortical white matter, not within the cortex. In other samples, we found suspected viral staining in locally ischemic regions (Fig. 8 D and Fig. S5, B andC). Remarkably, all the regions of positive viral staining showed no lymphocyte or leukocyte infiltration. This is in contrast to other neurotropic viruses (ZIKV, rabies virus, and herpes virus) in which the infection is typically accompanied by large number of immune cell infiltrates, including T cells. These findings suggest that, although SARS-CoV-2 has neurotropic properties and can infect neurons in patients, it did not invoke an immune response typical of other neurotropic virus in this particular case.

Discussion

We examined the potential for SARS-CoV-2 to infect neural tissues of both mice and human origin and demonstrate (B) DEGs from brain organoids infected with ZIKV [START_REF] Watanabe | Self-Organized Cerebral Organoids with Human-Specific Features Predict Effective Drugs to Combat Zika Virus Infection[END_REF] were compared with DEGs from SARS-CoV-2-infected organoids. (D) Enriched gene ontology terms (http://www.geneontology.org) for up-regulated genes from B. (E) Enriched gene ontology terms in SARS-CoV-2-infected cells (top) and SARS-CoV-2-negative bystander cells from 96-hpi organoids (bottom). (C) Heatmap of genes from metabolic pathways. (F) HIF1A staining of brain organoids that were mock infected versus 96 hpi. Scale bar = 100 µm. (G) Quantification of HIF1A-positive cells in SARS-CoV-2-infected organoids. Single-cell RNA-seq was performed in duplicate with one IPSC line (Y6) for reproducibility. HIF1A staining was performed with unique organoid, n = 4 per condition, from the same culturing batch, with images from n = 12 cortical regions with two IPSC lines, and Student's t test was performed (****, P < 0.0001). Experiments were performed twice for reproducibility. resuspended in PBS and aliquoted for storage at -80°C. Virus titer was determined by plaque assay using Vero E6 cells.

SARS-CoV-2 infection of organoids

Brain organoids in low-adhesion plates were infected with SARS-CoV-2 at multiplicity of infection 1.

SARS-CoV-2 infection (intranasal)

Mice were anesthetized using 30% vol/vol isoflurane diluted in propylene glycol. With a pipette, 50 µl of SARS-CoV-2 (3 × 10 7 PFU/ml) was delivered intranasally.

SARS-CoV-2 infection (intraventricular) Animals were anaesthetized using a mixture of ketamine (50 mg/kg) and xylazine (5 mg/kg) injected i.p. After sterilization of the scalp with alcohol and betadine, a midline scalp incision was made to expose the coronal and sagittal sutures, and a burr holes were drilled 1 mm lateral to the sagittal suture and 0.5 mm posterior to the bregma. A 10-µl Hamilton syringe loaded with virus and inserted into the burr hole at a depth of 2 mm from the surface of the brain and left to equilibrate for 1 min before infusion. Once the infusion was finished, the syringe was left in place for another minute before removal. Bone wax was used to fill the burr hole, and skin was stapled and cleaned. Following intramuscular administration of analgesic (meloxicam and buprenorphine, 1 mg/kg), animals were placed in a heated cage until full recovery. For the high condition, 5 µl of SARS-CoV-2 (3 × 10 7 PFU/ml) was used, and for the low condition, 5 µl of SARS-CoV-2 (3 × 10 6 PFU/ml).

Sample staining and iDISCO+ clearing

Whole-brain vasculature staining was performed following the iDISCO+ protocol previously described (Renier et al., 2016), with minimal modifications. All steps of the protocol were done at room temperature with gentle shaking unless otherwise specified. All buffers were supplemented with 0.01% sodium azide (Sigma-Aldrich) to prevent bacterial and fungal growth. Brains were dehydrated in an increasing series of methanol (Sigma-Aldrich) dilutions in water (washes of 1 h in methanol 20, 40, 60, 80, and 100%). An additional wash of 2 h in methanol 100% was done to remove residual water. Once dehydrated, samples were incubated overnight in a solution containing a 66% dichloromethane (Sigma-Aldrich) in methanol and washed twice in methanol 100% (4 h each wash). Samples were then bleached overnight at 4°C in methanol containing 5% hydrogen peroxide (Sigma-Aldrich). Rehydration was done by incubating the samples in methanol 60, 40, and 20% (1 h each wash). After methanol pretreatment, samples were washed in PBS twice, for 15 min and 1 h, in PBS containing 0.2% Triton X-100 (Sigma-Aldrich) and further permeabilized by 24-h incubation at 37°C in Permeabilization Solution (20% DMSO [Sigma-Aldrich] and 2.3% glycine [Sigma-Aldrich] in PBS-T).

To start the immunostaining, samples were first blocked with 0.2% gelatin (Sigma-Aldrich) in PBS-T for 24 h at 37°C, and the same blocking buffer was used to prepare antibody solutions. A combination of primary antibodies targeting different components of the vessel's walls were used to achieve continuous immunostaining. Antibodies to podocalyxin and CD31 were combined with antibodies against the nucleocapsid (N) from GeneTex. Primary antibodies were incubated for 10 d at 37°C with gentle shaking, then washed in PBS-T (twice for 1 h and then overnight) and finally newly incubated for 10 d with secondary antibodies. Secondary antibodies conjugated to Alexa Fluor 647 were used to detect podocalyxin and CD31, while the NUCLEOCAPSID protein was stained with a secondary antibody conjugated to Alexa Fluor 555. After immunostaining, the samples were washed in PBS-T (twice for 1 h and then overnight), dehydrated in a methanol/water increasing concentration series (20,40,60,80, and 100%, 1 h each, and then methanol 100% overnight), followed by a wash in 66% dichloromethane and 33% methanol for 3 h. Methanol was washed out with two final washes in dichloromethane 100% (15 min each), and finally the samples were cleared and stored in dibenzyl ether (Sigma-Aldrich) until light sheet imaging.

Light sheet imaging

We imaged with a 4× 0.35-NA objective cropped elongated field of view (600 × 2,200 µm) covering the narrow waist of the light sheet at 1.63 µm/pixel of lateral resolution, 1.6-µm spacing. A reference channel for registration to the annotated atlas using the sample autofluorescence was acquired at 5 µm/pixel. The acquisitions were done on a LaVision Ultramicroscope II equipped with infinity-corrected objectives. The microscope was installed on an active vibration filtration device, itself put on a marble compressed-air table. Imaging was done with the following filters: 595 nm/40 mW for Alexa Fluor 555, and 680 nm/30 mW for Alexa Fluor 647. The microscope was equipped with the following laser lines: OBIS-561 nm, 100 mW, OBIS-639 nm, 70 mW, and used the second-generation LaVision beam combiner. The images were acquired with an Andor complementary metal oxide semiconductor (CMOS) sNEO camera. Main acquisitions were done with the LVMI-Fluor 4×/0.3 WD6 LaVision Biotec objective. The brain was positioned in sagittal orientation, cortex side facing the light sheet. A field of view of 400 × 1,300 pixels was cropped at the center of the camera sensor. The light sheet numerical aperture was set to NA 0.1. Beam width was set to the maximum. Only the center-left light sheet was used. Laser powers were set to 100% (639 nm) or 10% (561 nm). Tile overlaps were set to 10%. The acquisition routine was set to the following: (1) Z-drive → save ome.tif stack; (2) filter change → Z-drive → save ome.tif stack; (3) change x position → repeat 1,2 12 times; and (4) change y position → repeat 1,2,3 6 times. At the end of the acquisition, the objective is changed to a MI PLAN 1.1×/0.1 for the reference scan at 488 nm excitation (tissue autofluorescence). The field of view is cropped to the size of the brain, the z-steps are set to 6 µm, and light sheet numerical aperture is set to 0.03 NA.

Viral RNA analysis

At indicated time points, mice were euthanized in 100% isoflurane. Brain tissue was placed in a bead homogenizer tube with 1 ml of PBS and 2% FBS. After homogenization, 250 µl of this mixture was placed in 750 µl Trizol LS (Invitrogen), and RNA was extracted with an RNeasy mini kit (Qiagen) per manufacturer protocol. To quantify SARS-CoV-2 RNA levels, real-time quantitative PCR (RT-qPCR) was performed using Luna Universal Probe Onestep RT-qPCR kits (New England Biolabs) with 1 µg of RNA and the Centers for Disease Control and Prevention real-time RT-PCR primer/probe sets for 2019-nCoV_N1.

Viral titer

Brain homogenates were centrifuged at 3,900 g for 10 min, and supernatant was taken for plaque assays. Supernatant at limiting dilutions was incubated on Vero E6 cells in MEMsupplemented NaHCO 3 , 4% FBS, and 0.6% Avicel RC-581. Plaques were resolved 48 hpi by fixing in 10% formaldehyde for 1 h followed by staining for 1 h in 0.5% crystal violet in 20% ethanol.

Stem cell culture

Human Y6 and Y1 iPS lines were obtained from Yale Stem Cell Center and New York Cell Stem Foundation, respectively. Cells were verified as being pluripotent, having normal karyotype, mycoplasma free, and cultured in feeder-free conditions on Matrigel-coated plates with mTeSR Plus culture medium (Stem Cell Technologies) and passaged using ReLeSR (Stem Cell Technologies).

Teratoma formation 10 6 Y6 iPSC cells were collected by collagenase treatment and resuspended in 100 ml of DMEM/F12, collagen, and Matrigel mix (2:1:1 ratio). Cells were injected i.m. into immunodeficient Rag2 -/-GammaC -/-mice. After 8 wk, teratomas were harvested, fixed, and subjected to paraffin embedding and H&E staining.

G-band staining for karyotype analysis Y6 iPSC small clumps were seeded on glass slides precoated with Matrigel and fed with mTESR medium for 3 d. Media were then switched to DMEM basal medium supplemented with 10% FBS for another 3 d, and the slides were transferred to Yale Cytogenetic Laboratory for G-band staining.

Neural progenitor cell culture Y6 iPS lines were differentiated to NPCs on Matrigel-coated plates using the monolayer protocol of the StemDiff SMADI Neural Induction kit (#08581; Stem Cell Technologies) for two passages and then maintained in StemDiff Neural Progenitor media (#05833; Stem Cell Technologies). 12-d NPCs were used for all experiments.

Cerebral organoid culture

For preparation of embryoid bodies, 9,000 single cells were seeded in each well of low-attachment 96-well U-bottom plates and 10 µM Y-27632 ROCK inhibitor for 1 d. Cerebral organoids were generated following exactly the previously established protocol [START_REF] Lancaster | Generation of cerebral organoids from human pluripotent stem cells[END_REF][START_REF] Lancaster | Cerebral organoids model human brain development and microcephaly[END_REF], using an orbital shaker for agitation.

Immunostaining

Brains of infected mice and organoids were collected and fixed in 4% paraformaldehyde (PFA). Samples were then dehydrated in a 30% sucrose solution. Cryostat sections were blocked in 0.1 M Tris-HCl buffer with 0.3% Triton and 1% FBS before staining. Slides were stained for IBA-1 (nb100-1028; Novus Biologicals), GFAP (ab4674, [Abcam] and rabbit anti-SARS-CoV-2 nucleocapsid [GeneTex]). Slides were mounted with Prolong Gold Antifade reagent (Thermo Fisher Scientific). All slides were analyzed by fluorescence microscopy (BX51; Olympus).

Organoids were fixed in 4% PFA in a BSL3 facility then moved to 30% sucrose solution at 4°C for ≥24 h. Organoids were embedded in tissue freezing medium, cut into 20-µm sections using a cryostat, and mounted on slides. After blocking and permeabilization with 0.25% Triton X-100 and 4% donkey serum, sections were incubated overnight with primary antibody: ms anti-Pax6 (1:300, #561462; BD PharMingen), rb anti-ACE2 (1: 500, ab15348; Abcam), ms anti-Sox2 (1:500, sc-365823; Santa Cruz), rat anti-CTIP2 (1:500, ab18465; Abcam), rb anti-TBR1 (1: 500, ab31940; Abcam), ms anti-MAP2 (1:500, MAB3418; Millipore), rb anti-SARS-CoV-2 nucleocapsid (1:250, GTX635679; GeneTex), and rb anti-HIF1 α antibody (GTX127309; Genetex). TUNEL assay was performed using the Click-It Plus TUNEL Assay for In Situ Apoptosis Detection (C10617; Thermo Fisher Scientific) per the manufacturer's instructions. Alexa Fluor 488, 555, and 647 antibodies were applied for 1 h at room temperature (1:500) after three 10-min PBS washes. To mark nuclei, DAPI (1: 3,000) was added to the secondary antibody incubation. Slides were then washed three times in PBS and then mounted with VectaShield Anti-Fade Mounting Medium. Images were acquired using an LSM 880 confocal microscope (Carl Zeiss) and prepared using Fiji (National Institutes of Health).

NPCs on coverslips were fixed in 4% PFA for 15 min at room temperature in BSL3 conditions and washed three times with PBS. After blocking and permeabilization with 0.25% Triton X-100 and 10% donkey serum, coverslips were incubated overnight with primary antibody as described above. TUNEL assay was performed as described above. Alexa Fluor 488, 555, and 647 antibodies were applied for 1 h at room temperature (1:500) after three 10-min PBS washes. Coverslips were then washed three times in PBS and then mounted on slides with Prolong Diamond Anti-Fade Mountant with DAPI (Thermo Fisher Scientific). Images were acquired using an LSM 880 confocal microscope (Carl Zeiss) and prepared using Fiji.

Human FFPE sections were heated for 30 min at 60°C and treated with xylenes followed by rehydration in decreasing concentrations of ethanol (100, 90, 80, and 70%). Antigen retrieval was performed using a pressure cooker (Biocare Decloaking Chamber) by boiling in sodium citrate (pH 6.0) for 20 min at 115°C. After blocking and permeabilization with 0.25% Triton X-100 and 4% donkey serum, sections were incubated overnight with primary antibody: rb anti-ACE2 (1:500, ab15348; Abcam) and ms anti-NeuN (1:500, ab104224; Abcam). Alexa Fluor 488 and 555 antibodies were applied for 1 h at room temperature (1:500) after three 10-min PBS washes. To mark nuclei, DAPI (1:3,000) was added to the secondary antibody incubation. Slides were then washed three times in PBS to remove detergent and treated with TrueBlack Lipofuscin Autofluorescence Quencher (Biotium) according to the manufacturer's instructions. Sections were mounted with

General Conclusion

The vasculature, an intricate and expansive network, plays a pivotal role in the brain, interweaving with other neural networks and cell populations. Its primary mission is to ensure the efficient distribution of vital nutrients and oxygen, essential for energy production. Achieving this delicate equilibrium involves a multitude of mechanisms, known as neurovascular coupling interactions. Disruptions in this delicate balance are suspected culprits in many cerebral vascular diseases, including neurodegenerative conditions and neuropathies. They also feature prominently in the regulatory pathways governing neural apoptosis during events like strokes.

The precise mechanisms governing vascular growth, rearrangement, neurovascular interactions, and related phenomena remain enigmatic, particularly in terms of their scaling effects, topological characteristics, and underlying molecular processes. To probe these phenomena and to understand the integral role of vasculature within the broader cognitive system, it is imperative to acquire large-scale vessel network datasets alongside molecular and functional mappings. Regrettably, this represents a substantial technical challenge. In response to this challenge, we have developed a comprehensive pipeline capable of facilitating whole-organ vascular imaging within a reasonable timeframe. Additionally, we offer a versatile analysis toolkit designed to characterize the state and evolution of the vascular system under various pathological conditions, including strokes and hearing loss. Moreover, we propose to monitor the developmental trajectory of the vascular network in healthy subjects and under various pathological conditions. To further support the research community, we have made available resources in the form of atlases and computational toolkits that facilitate the study of pup development.

This comprehensive approach ultimately leads to the full characterization of the vascular system, revealing substantial disparities in terms of density and topology. Particularly, I found that in cortical regions, these disparities are sufficiently significant to infer the functional attributes of different cortical areas, whether integrative or somatosensory in nature. Additionally, we have explored the plasticity of vasculature in pathological conditions, such as the rewiring patterns and time-lapse observations following strokes, as well as vascular degradation in congenitally deaf mice. Our findings provide compelling evidence of the correlations between neural activity, stress, and vascular topology.

80

In further pursuit of understanding vascular heterogeneities, I have studied the evolution of vascular characteristics throughout development. Vascular development unfolds in distinct waves, manifesting at different time intervals across meta regions. Notably, somatosensory systems, both cortical and thalamic regions, appear to co-develop, suggesting a link between the length of the developmental timeframe and the observed orientation/density heterogeneities.

Early perturbations in the system, such as sensory deprivation, result in alterations in adult vascular topology. This implies that reduced neuronal activity during development leads to lower vascular density in the affected regions. These findings suggest the co-existence of activity-dependent and pre-patterning signaling pathways controlling vascular development. It would be important to follow up on this work to better characterize vascular plasticity in the adult brain.

Figure 1 :

 1 Figure 1: Examples of natural and artificial networks

Figure 2 :

 2 Figure 2: brain size across mammalian species, from Herculano-Houzel et al 2020. Among all the above species, brain sizes and neuronal densities vary by several orders of magnitude. The cellular density of the cerebral tissue is most likely an important factor dictating vascular organization and density.

Figure 3 :

 3 Figure 3: Organization of the neuro-vascular unit -adapted from Zhao et al 2015. Neuronal and vascular networks are tightly coupled through the control of hemodynamic and nutrient supply. They interact via support cells, such as the pericytes and astrocytes, controlling the exchange of solutes, nutrients, and blood flow.

Figure 4 :

 4 Figure 4: Li, W. et al. Cerebral vascular alterations in a mosue model of type 2 diabetes. Changes in the composition of the blood affects directly the vascular topology and could createweaknesses in the network.(Diabetes 59, 228-235 (2010)).

Figure 5 :

 5 Figure 5: effect of sars-cov infection on mouse brain vasculature, adapted from Wenzel et al 2021. This viral infection is believed to cause the appearance of "string vessels" that are not functional. Sars-Cov 2 has been linked to brain alteration in both mice and humans in some cases, assessing the importance of the vascular component in brain diseases.

Figure 6 :

 6 Figure 6: de Vinci cardiovascular circuit -1508. The description of the vascular system as a closed circulation network has revolutionized our understanding of general physiology.

Figure 7 :

 7 Figure 7: adapted from Risau et al 1997. Initial formation of the Vascular network during early stages of development. Precursor mesodermic cells, once differentiated into vascular cells, aggregate together to form "blood islands", that will under the influence of various growth factor, bind together until they form a primitive capillary plexus. Later on during the development this primitive network will rewire and complexify via sprouting and pruning processes influenced by molecular gradients.

Figure 8 :

 8 Figure 8: from Thion, Garel et al 2017. Several waves of developmental processes occur in the embryonic and post-natal brain. Interestingly, most of the cerebral angiogenesis occurs in the mouse in the first week of postnatal life. This is at the same time as neuronal activity shapes neuronal connectivity via synaptogenesis and synaptic pruning, suggesting a strong interaction between the refinement of neuronal and vascular networks.

  To a d d r e s s t h e d e f o r m a t i o n s v i s i b l e o n b r a i n s e c t i o n s , b l o c k f a c e i m a g i n g h a s b e e n i n t r o d u c e d .

  The principle of light sheet microscopy is quite simple: the transparent tissue is put in a container filled with a solution of the same refraction index as the sample. The sample is then illuminated from the side by a laser forming a thin planar Gaussian beam. Illuminated stained tissues emit fluorescent light that is then collected by an orthogonal microscope objective lens located above the sample (Fig 9). Depending on the fluorophores, different excitation wavelengths can be used. For very large objects open-top light sheet designs use a setup tilted to 45 degrees above the samples (usually a large slice of primate tissues). Since the gaussian beam has a narrow field allowing isotropic and in-focus imaging, data are acquired as thin tiles centered on the waist of the light sheet to get an even resolution. Columns of data are thus acquired and later stitched together using stitching pipelines.

Figure 9 :

 9 Figure 9: Dodt 2007: Principle of ultramicroscopy. (a) The sample is illuminated from the sides by a laser forming a thin sheet of light from low numerical aperture cylindrical lenses. Light thus only illuminates a thin optical section and fluorescence emission is collected by a perpendicular high numerical aperture objective lens. The image is projected through the tube lens onto the camera target. (b) Photograph of a whole mouse brain in an organic solvent clearing solution. The highly autofluorescent myelin is visible.

Figure 10 :

 10 Figure 10: from Seo et al 2016 Tissue clearing techniques. (A) Light scattering in biological tissues can be reduced by removal of lipid and Refraction Index (RI) matching. (B) Simple immersion in a high-RI aqueous solution renders the tissue modestly transparent by homogenizing scattering throughout the sample. (C) Delipidation and dehydration/hyper-hydration followed by refractive index matching. (Top) For solvent-based clearing, the tissue is incubated in dehydrating solvent for delipidation and is moved to a high-RI clearing solvent where RI matching and additional delipidation occur. (Bottom) The sample is placed in an aqueous solution that contains high concentration of non-ionic detergent and denaturant, where delipidation, hyperhydration, and RI matching take place. (D) A biological sample is first transformed into a tissue-gel hybrid by hydrogel embedding (Top) or glutaraldehyde fixation (Bottom), where the gel network increases the tissue integrity. The tissue-gel hybrid then can withstand extensive delipidation by incubation in ionic detergent (SDS) assisted by electrophoresis or heating.

Finally, organic

  solvent-based methods like 3DISCO (Erturk et al., 2012) are cheap and have very high transparencies. Their drawback is the toxicity and the difficulty of preserving the fluorescence of genetically encoded fluorophores such as GFP or tdTomato.

Figure 11 :

 11 Figure 11: adapted from Arganda-Carreras et al 2017. This caption illustrates the generic image segmentation pipeline strategy, which Tubemap follows. Any segmentaion requires pre and post processing of the data. These can vary depending of the specificities of the image acquisition system.Most preprocessing steps aim to enhemce the signal on noise ratio, remove bckground noise, increase contrast and generally get rid of any system artifacts introduce during acquisition. The quality of the pre processing will stronglu affect the final results after segmentation. The segmentaion part itself can involve many different strategies. One main strategy consists in enhencing specific trait of the object to segment (contour etc) followed by a thresholding. Another common strategy is to use iterative algorithm like watershed, that will iteratively increse the segmented region following some cost function set by the user. Finally postprocessing aim to clean up the final resulting data. Sometime it can be manual curation of false positive, or simply filling up and smoothing the segmentation outputs. The algorithm use once again strongly depends of the type of object to be segmented. In this regard TubeMap is quite a generic pipeline, and can possibly be adapted to many different type of data.

Contour detection :

 detection It allows the estimation of the boundaries of an object based on the image's first or secondorder derivative. It usually requires additional processing to obtain continuous borders and can typically be combined with watershed segmentation algorithms. Watershed transform: It uses an analogy between grayscale images and landscape topology, where highintensity pixels are peaks and low-intensity ones are valleys. It then mimics water flooding through this landscape and builds dams when waters from different sources meet. Each separate pool of water defines an object. It works particularly well with objects that have bright edges. Pixel classification: Binarization can be seen as a pixel classification task where pixels must be classified together with other pixels presenting similar intensity or localization. Classification requires the definition of a metric function. The most popular unsupervised method is k-means clustering. There are also supervised methods where the user provides training data that is already labeled.

Figure 12 :

 12 Figure 12: Architectures for image reconstruction, adapted from Belthangady 2019 The U-Net is a variation of the encoder-decoder network. It adds shortcuts or skip connections between the encoder and decoder branches. These additional connections are beneficial for image restoration tasks because they send fine image details directly to the decoder path, skipping the bottleneck of the deep layers. For image restoration tasks where the identity function is a good first guess, residual connections that add the input image to the output image can be used to learn residual mapping. Overall, skip and residual connections have been shown to help train deeper networks by preventing vanishing gradients. (Olaf Ronneberger, Philipp Fischer, 2021)

  To m a k e i t s c a l a b l e f o r l a r g e p h y s i c a l o r b i o l o g i c a l n e t w o r k s , v a r i o u s m e t h o d s h a v e b e e n u s e d . M a r k o v C h a i n Monte Carlo (MCMC) was first proposed since it is straightforward to implement (Peixoto, 2014) . The MCMC algorithm uses sampling to evaluate the expected value of a random variable (Fig 14).

Figure 14 :

 14 Figure 14: from Peixoto et al 2014. This table summarizes the main variations of the Stochastic Block Model algorithm. Those variants may consider the degree distribution of vertices in each cluster to avoid clustering high degree vertices which is a shortcoming of the standard SBM method. In other variants, the user can consider edge weigh to control the clustering. Another interesting feature is the possibility to set the number of clusters. This is particularly important to benchmark for modularity, or if many modules are needed. Hierarchical SBM perform standard SBM at different scales recursively. All these techniques have their pros and cons. In the case of this work, the standard SBM was used to avoid biasing the module search.

Equation 4 :

 4 Poiseuille's Law and formulation of the flow

  Article 1: Mapping the fine scale organization and plasticity of the brain vasculature Mapping the fine scale organization and plasticity of the brain vasculature, Cell 2020 Kirst C, Skriabine S, Vieites-Prado A, Topilko T, Bertin P, Gerschenfeld G, Verny F, Topilko P, Michalski N, Te s s i e r -Lavigne M, Renier N. Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature.Cell. 2020 Feb 20;180(4):780-795.e25. doi: 10.1016/j.cell.2020.01.028. Epub 2020 Feb 13. PMID: 32059781

. 1

 1 Figure 1. Deep Immunolabeling and Imaging of the Vasculature with Light Sheet Microscopy (A-E) Complementary immunolabeling against CD31, podocalyxin, and Acta2 enables deep immunolabeling of the vasculature in intact perfused brains. (A) Immunolabeling strategy.

  (B) Light sheet scanning strategy along the narrow waist to accelerate acquisitions and maximize the optical resolution. (C) Whole brain immunolabeling of the vasculature with iDISCO+ for CD31 and podocalyxin (blue) and Acta2 (red), imaged with light sheet microscopy. Transverse cortical 50-mm thick projection plane from a central tile. (D) Acta2 enables the discrimination of arteries in the brain parenchyma: an Acta2+ artery is visible (a), and a nearby Acta2-vein (v). (E) Summary of the signal distribution in the brain. (F) Light sheet scans produce non-stitchable columns. A projection along the overlapping region of the tile is shown with misalignments noticeable in a few planes at the center of the inset. (G) WobblyStitcher addresses recurrent stitching issues to produce seamless datasets. The same alignment is shown after non-rigid alignment of the planes, correcting the vessel duplications. (H) Presentation of the stitching pipeline: (1) tile pairs are first aligned in z; (2) each plane is aligned individually along x,y; (3) tracking of the alignment quality for a smooth correction of the rigid stitching artifacts and smoothing of the resulting trajectories; (4) global optimization of all plane positions. (I) Raw data resulting from the labeling, imaging, and stitching strategy, shown in ventral, dorsal, and lateral orientations. A thick projection of part of the middle cerebral artery territory is shown on the right, with striatal penetrating arteries visible. Scale bars represent 300 mm in (C), (D), (F), and (G) and 1 mm in (I).

Figure 2 .

 2 Figure 2. Vasculature Segmentation through Multi-path Binarization, Tube Filling, and Parallel Skeletonization (A) Presentation of the image-processing pipeline for the binarization of vascular immunolabeling. Original data are on the left, and the binarization result is on the right. (B) The multi-path TubeMap binarization is adaptable to segment complex immunolabeling signals and contains steps optimized for large 3D datasets (green circles) and steps specifically designed for the vasculature (yellow circles). (C) A novel strategy for fast removal of light sheet striping artifact using a combination of 3D-rank filters. The input shows stripe artifacts (left). A local percentile filter with elongated structural element (red box) estimates the light sheet artifact at each voxel (middle left). To avoid removing vessels in that orientation it is compared to an estimate of the overall vessel intensity using a cubic structural element (brown box, middle right). The result is subtracted to give the corrected image (right). (D) Pseudo deconvolution method to correct for shadows diffusing the edge of large vessels (top left). High threshold pixels are determined (top right) and a blurred version subtracted from the original image to give a pseudo-deconvolved image (bottom left), the thresholded version (bottom right) separates the large bright branches. (E) Result of the 3D equalization (left) of the light sheet corrected image in (C). The vessel intensities are homogenized. (F) Result of the Hessian tube filtering to enhance smaller and weaker vessels. (G) Sum of the different multi-path filters. (H) Convolutional neural network applied to the binary mask containing hollow vessels: a single plane is presented, the 3D projection of the region, and the resulting graph. The network fills empty tubes recognized from the binary mask. (I) Fast centerline extraction via parallel thinning.. The centerline and binarized vasculature are then used to build a graph of the vasculature network and the geometry of the vessels.

Figure 3 .

 3 Figure 3. Unsupervised Construction of Labeled Vascular Datasets (A) The full TubeMap pipeline.(B-F) Slices through a vascular brain graph, obtained from a C57bl6 3-month-old female mouse, plotted with ClearMap 2.0. Arteries are indicated in red, putative veins in blue, and capillaries are colored according to their regional annotation. (B) 250-mm coronal and sagittal slices of the graph, showing the general organization of the cerebral vasculature and its regional variations. The inset shows a detail at the level of the hippocampus, with an arterial arc visible. (C) 750-mm coronal and sagittal slices taken at the same levels, showing only labeled veins and arteries. Insets show sagittal and coronal views of the hippocampus, showing alternating arterial and venous arcs. (D) 250-mm projection at the level of the brainstem midline, showing the radial organization of veins and arteries in the hindbrain and superior colliculus. (E) Detail of a 1-mm midline sagittal slice at the level of the cerebellum: veins are predominantly located in the white matter and arteries at the pial surface of the lobules. (F) Detail of the anterior commissure at a midline level (sagittal section). The white annotated vessels (fiber tracts) are longer and less dense that the gray matter vessels (red, thalamus; light blue, striatum). Abbreviations: SO, stratum oriens, SP, stratum pyramidale, SLM, stratum lacunosum moleculare, SM, stratum moleculare, SG, stratum granulosum, H,: hylus. Scale bars represent 1 mm, except in (F) (250 mm).

Figure 4 .

 4 Figure 4. Properties of the Whole Brain Vascular Graph (A) Added vessel lengths in meters across major brain regions in 1 hemisphere (top, example of 1 brain), and average length density across all brain regions in meters per cubic mm (bottom, n = 3). The distribution of lengths densities is highly variable across brain regions. A detail of cortical layers is shown as an inset. (B) Branch-point density across brain regions, showing similar trends as in (A). (C) Spatial visualization of the distribution of branch point densities across the brain revealing hotspots of high vascular density related to specific nuclei. (D) Vessel density versus branch point density are highly correlated across brain regions (r 2 = 0.94). The proportion of vessels to branch points of 1.78 indicates a systematic redundancy in the network. (E) Counts of 4-point motifs in the vasculature graph throughout the brain. Although the fork motif dominates, motifs with local cycles constitute $15% of the network. (F) Distribution of vessel branch lengths in major brain regions. (G) Joint distributions of branch length and radius across brain regions. (H) Example of a shortest path (in color) between an artery (red) and a vein (blue) with other possible paths (in gray). (I) Number of branches from all capillaries to their nearest artery or vein. (J) Length of the shortest paths from each capillary branch to the nearest artery (red) or vein (blue) and from arterial vessels to the nearest vein (light blue) and venous vessels to the closest artery (light red). (K) Joint distribution of the shortest path lengths from each capillary to their nearest artery and vein for major brain regions. (L) Roots of individual penetrating arteries isolated in the cortex. (M) Distribution of the sizes of arterial domains defined by shortest distance to the arterial roots in the cortex. Examples of small, medium-sized, and large domains are shown. Scale bars represent 1 mm.

Figure 5 .

 5 Figure 5. Regional Specificity of the Layered Structure of the Cortical Vasculature (A-C) Features extraction of the vasculature from cortical regions. (A) Branch point density. Higher densities are seen in the sensory cortices, in particular in layers 4. (B) Arterial branch points with two representative regions are shown: barrels (SSp-bfd) and primary motor (MOp). (C) Color-coded branch orientations show layering in the somatosensory and motor cortices. The vessels in the somatosensory cortex have predominant radial orientations in upper layers 2/3 (blue), while the motor cortex has a mix of tangential (green-pink) and radial orientations across all layers. (D) Distribution of features across layers (examples of the SSp-bfd and MOp, n = 3 brains).(E) Distances between the feature distributions (earth-movers distance between the normalized feature distributions). The clustering of the distances reveals 2 major groups of regions, containing preferentially (1) motor and integrative regions (red), and (2) primary and secondary sensory regions (green). Scale bars represent 1 mm (A), 500 mm (B) and (C), and 200 mm (C, insets). See TableS1for abbreviations.

Figure 6 .

 6 Figure 6. Application of TubeMap to the Study of Ischemic Strokes (A) Projection of the Acta2 labeling in a control brain or brains recovering for 3 weeks after subjected to an electrocauterization of the middle cerebral artery. Sagittal and ventral projections are shown. The insets show the destruction of the branch (yellow circle) connecting the middle cerebral artery (MCA) dorsal arbor to the circle of Willis. (B and C) Reconstruction of the vasculature in control (B) and stroked (C) brains. A 250-mm thick sagittal slice through the graph are shown with no noticeable difference in the position and numbers of major veins (in blue) or arteries (red). (D) Inset from the sliced graph at the level of the supplemental somatosensory (SSs) and mouth primary somatosensory (SSp-m) cortices. Reorientations of the capillaries across layers toward the center of infarct (dotted line) are noticeable in the stroked brain. Color-coding of the vessel orientations reveal a mix of radial (green) and planar (blue and red) orientations in controls, while planar (along the cortical plane) orientations are predominant in the stroked cortex. (legend continued on next page) Cell 180, 1-16, February 20, 2020 11

Figure 7 .

 7 Figure 7. Brain-wide Reorganization of the Vasculature in a Model of Deafness (A) Sliced graph in coronal orientation at the level of the primary auditory cortex in controls (left) and otoferlin mutant mice (right). Graphs of the auditory cortex arteries are shown. A reduction in the vessel density in the layer 4 of the auditory cortex in the mutant is visible in the mean voxelized branch densities (p value map: p < 0.01, Mann-Whitney test, n = 3). (B) Sliced graph at the level of the somatosensory barrel cortex. A higher density of vessels is visible in the mutant barrel cortex in the mean voxelized branch densities (n = 3, p < 0.01, Mann-Whitney test). (C) Mean branch densities and graphs of subcortical sensory relays are shown for the inferior colliculus, superior olivary complex, nucleus of the lateral lemniscus, and cochlear nucleus. (D) While the cochlear nucleus shows a reduction in branch densities, the somatosensory trigeminal principal nucleus (PrV) is not significantly different between control and otoferlin mutant mice. (E) Sub-domains of vessels inferred in the auditory cortex (AUDp) or barrel cortex (SSp-bfd) of control mice using a stochastic block model (SBM). A sub-domain represents a group of vessels stronger connected to each other than to vessels from another sub-domain. Vessels belonging to the same domain are colored in matching colors. (F) Views of the domains in layers 5 and 4 of the AUDp and SSp-bdf, respectively, in control or mutant mice. (G) Number of SBM sub-domains and vessel branches in mutant and control mice for AUDp and SSp-bdf. Results are presented across all cortical layers, or only layer 5 (L5) or 4 (L4). Otoferlin mutant mice show the same number of domains on average than controls, but the number of edges per domain changes in different regions (n = 3, Mann-Whitney test). Scale bars are 500 mm.

  anti-TRANSGELIN/SM22 Dilution 1:1500 Abcam Cat# ab14106, RRID: AB_443021, Lot: GR3247382-1 Rabbit anti-TRANSGELIN/SM22 Dilution 1:300 Proteintech Cat# 10493-1-AP, RRID: AB_219936, Lot: 00040337 Rabbit anti-CollagenIV Dilution 1:1000 Abcam Cat#ab19808 RRID: AB_445160, Lot: GR261208-23 Rabbit anti-Human Von Willebrand Factor Dilution 1

  ðI i ðxÞÀI j ðx + sÞÞ 2 where N ij ðsÞ =

Figure S7 .

 S7 Figure S7. Reconstruction of the Spinal Cord Vascular Network, Related to Figure 4

  itk snap Based on the paper atlas of P0, and P6 (Atlas of the developing mouse brain by Halasz et Tsalis, 2007) and Young et al 2020 paper, we annotated by hand 57 different brain regions using itk-snap open source software. The list of the annotated region is the following: AUDp, AUDs, SSp, PAG, PrV, VllN, Striatum, VISp, OLF, AmyC, AmyS, SpV, RT, A, IntC, HYP, MB, septums, Fasciculus Retroflexus, Stria Terminalis, Fornix, Bed nucleus of the stria terminalis, Fimbria,ventral hippocampal commissure,Cortical subplate, Globus pallidus,Neocortex, HIP, TH, CB, IC, SC, MOp, MOs, SSp-bfd, VPM, VPL, MNTB, SOC, VN, IIIN, IVN, VIN, Deep cortical tracts,Lateral ventricle, MGN, vLGN, Substantia nigra reticular, Red Nucleus, medial habenula, lateral habenula, Antéro Dorsal nucleus of the Thalamus, AnteroVentral nucleus of the Thalamus ,Inferior Olive.

  More formally: validation We ask experts in anatomy to check our segmentation and compare with the current resources available (Allen brain atlas of adult mice, Halasz et Tsalis, 2007, and Young et al 2020. We also created an open source free to browse resource repository where users are invited to explore the data and propose modifications based on their knowledge, needs in terms of annotated regions. l atlas labelling the output of the previous part are a collection of autofluorescence templates and their corresponding annotation files at the resolution 25 um × 25 um × 25 um for the timepoints P1, P3, P5, P7, P9, P12, P14, and P21.

  Figure 1: A the atlas alignment and averaging pipeline. 20 samples were used to produced each of the atlas at every timepoint B The manual annotation pipeline, and 3D laplacian smoothing strategy. the annotation was done on the coronal view

  As a test, I worked on training a novel DNN specifically on very large vessels (Sup Fig 1) to improve the detection of those large surface veins and arteries. While this network improved the detection of larger vessels, further characterization is needed to fully validate the use of this new model.

  figure 17: (A) pipeline summary (adapted from Kirst et al 2020) (B) antibodies conjugation strategy: the diffusion of the antibodies through the tissues is a limiting factor that results in high heterogeneities in the quality of the tissue labeling, with higher signal values on the surface than on the brain center, with a decreasing signal to noise ratio. Conjugating the primary antibodies allows a better diffusion in the tissue than the use of a primary-secondary assembly, and increases signal to noise ratio in deeper regions, leading to better data quality. (C) current vessel filling DNN results as described in (Kirst et al 2020). This DNN was specifically tuned to detect most medium to low size empty vessels and fill them up to allow the 3D thinning step of the TubeMap pipeline to run correctly. However very large vessels were not corrected. First because of their size they would not fit into the training 3D data chunk. Increasing the data chunk during training require modern GPU with large amount of GPU RAM, and slows down the training. Second, only a few vessels have large sizes, representing a very small proportion of the total vessels counts, making it harder to tune a neural network to recognize them. I trained another DNN that specially target large vessels for filling. It runs in addition to the previous DNN and relies on data downsampling to detect those large vessels. (D) new DNN for very large vessels filling result and comparison. [prior to DNN filling processing neither the basilar artery or the smaller arterioles and veinules are filled. After the first DNN, medium sized arterioles and veinules are filled. After the second DNN even the Basilar artery is filled.

  figure 18: (A) Stochastic Block Model principle illustration. The stochastic block model subdivides graph into communities, i.e subsets of nodes with higher connective density with one another than the rest of the network. It is an important statistical model in network science for discovering community structures in complex graphs.

  figure 19: (A-B) The auditory cortex is a region composed anatomically of several subregions (primary, ventral and dorsal), themselves subdivided in layers. The neuronal substrate of these different part is involved in different processing circuitry. Several studies link community structure and role, which is why SBM on region like AUD is interesting to unveil. The SBM here seem to reveal some upper-and lower-layers community structures. Specifically, it is useful as a benchmark for the quality of a network. An altered network is likely to present dissimilar community property, such as an altered modularity or a significantly different number of vertex. In this panel we compare the community partitioning of otof-/-congenitally deaf mice auditory cortex with the auditory cortex of control mice. The lower left plot show the vascular density and orientation profiles of Otof-/-and control mice at 2 months and 6 month of age. The y axis represents respectively the proportion of radially oriented vessels on the left and number of vessel branch point on the right. The x axis represents the cortical depth, from the surface (at 0) to the last cortical layer. Starting from layer 4 down to the lower cortical layers, otof-/-mice show a total of 12 % less vessels than controls in auditory areas at 2 months, and this difference increases even further with aging (6 months), even though vascular pruning is a normal mechanism of vascular aging. The lower right panel shows the statistical difference in the modularity and number of modules between controls and otof-/-mice. Both are statistically lower in otof-/-mice. These modifications are mostly driven by the decrease of the vascular density in otof-/-in the auditory cortex. Finally we show what the auditory cortex and partitioning would look like after removal of 40% random vessels. (C) In mathematical graphs, removing an edge whose vertices are of degree 3, leads to topological degree 2 vertices that are biologically irrelevant. Therefor an algorithm allowing edge pruning without affecting the degree distribution of the graph vertices is necessary to model the vessels depletion the vasculature undergoes in otof-/-auditory cortex. This panel presents a description of the pruning algorithm we developed to mimic vessels removal. To quickly summarize, remaining degree two vertices are remove and replace with a single edge at each prunig iteration. Vessel removal is considered as random but probability weight are put on edges based on different criteria, such as edge length, radius, extracted fraction, and distance to surface. (D) Evolution of a vascular graph modularity and number of modules when edges are progressively removed. Healthy graphs reach pathological values when 12 % of the vessels are removed which corresponds to the average differences we witness in hearing loss auditory cortices and control auditory cortices.

  To t e s t t h e e ff e c t o f n e u r o n a l a c t i v i t y o n t h e f i n a l v a s c u l a r o r g a n i z a t i o n d u r i n g d e v e l o p m e n t , I u s e d a m o d e l o f congenital deafness. The Otoferlin mouse line presents a mutation of the protein responsible for synaptic

( C )

 C figure 22: (A) Orientation and BP profiles of the nose and barrels fields somatosensory cortical area in control and whisker deprived mice. The y axis represents respectively the proportion of radially oriented vessels on the left and

  is a tool based on cloning, allowing different developers/users to work in parallel on the same document, without altering the work of the others. The chronology of the document is saved which easily allows to retrieve back a previous version of the work, preventing the loss of an already saved work. It is widely used in programming projects typically when several people are working on the different facets of the same project.Streamlining the sample preparationAs previously described in the paragraph about the comparison between the different clearing techniques, the iDISCO protocol used in this thesis, is highly toxic. It involves manipulating methanol, DCM, and DBE which are all highly corrosive and potentially cancerogenic for humans. It is also quite tedious since it requires changing the sample solution every hour, making it difficult to work in parallel on another experiment. A last limitation resides in the step-like nature of the hydration/rehydration steps of the iDISCO protocol. The idea is to slowly shift from an aqueous to a 100% methanol medium, requiring a gradual change of concentration. Handily it means increasing the concentration of one of the solvents by 20 % every hour. Instead of a smooth gradient shift, the concentration profile is rather step-like. It makes the full protocol

  , top panel). Conversely, SARS-CoV-2-negative cells showed a mitochondrial catabolic state with the up-regulation of alcohol metabolism, cholesterol synthesis, and regulation of cell death (Fig.5 E, bottom panel). These two cell populations displayed antagonizing pathway enrichment, with the infected cells responding to hyperoxia and the bystander cells responding to hypoxia (Fig.5 E, highlight in yellow). The hypermetabolic state is unique to the SARS-CoV-2-infected cells (Fig.5 C)a n d highlights the ability of SARS-CoV-2 to hijack the host neuron machinery to replicate (Fig.2 D). Finally, we confirmed that infection by SARS-CoV-2 induced a locally hypoxic environment in neuronal regions by staining for HIF1α (Fig.5, F and G)i n mock-infected and SARS-CoV-2-infected organoids. Together, these results indicate the potential of SARS-CoV-2 in manipulating host metabolic programming, which may create a resource-restricted environment for cells.

Figure 1 .

 1 Figure 1. Pluripotency of Y6 line. (A) Sample images of cell types of three germ layers in teratomas following transplantation to Rag2 -/-GammaC -/-mice. Scale bar = 100 µm. (B) Sample image showing normal karyotype for Y6 line. Scale bars = 100 µm.

Figure 2 .

 2 Figure 2. SARS-CoV-2 infects human brain organoids and induces cell death. Human brain organoids were infected with SARS-CoV-2 and collected 24 or 96 hpi to analyze for different cellular markers. (A) Images of brain organoids looking at SARS-CoV-2 infection (in red) 24 hpi (see Fig. S2 C for additional images). Scale bar = 200 µm for zoomed-out images and 25 µm for zoomed-in images. (B) Images of brain organoids looking at SARS-CoV-2 infection (in red) 96 hpi (see Fig. S2 D for additional images). Scale bar = 200 µm for zoomed-out images and 50 µm for zoomed-in images. (C) Quantification of SARS-CoV-2-positive cells in a single microscope image of cortical region of organoids (A and B). (D) Tiled image of 96-hpi organoid. (E) Electron microscopy image of SARS-CoV-2 viral particles in brain organoids (see Fig. S3 for uncropped and additional images). Scale bar = 100 nm. (F) Organoids were stained with TUNEL to evaluate cell death at 96 hpi. Scale bar = 100 µm. (G) Quantification of SARS-CoV-2 and TUNEL double-positive (yellow) or SARS-CoV-2-negative, TUNELpositive (green) cells over total TUNEL-positive cells. (H) Quantification of SARS-CoV-2 and TUNEL double-positive (yellow) or SARS-CoV-2-positive, TUNELnegative (red) cells over total SARS-CoV-2-positive cells. (I) Correlation between the frequency of TUNEL-positive cells and presence SARS-CoV-2 in different regions of the organoid. (J) Representative image of TUNEL and SARS-CoV-2 staining showing high-density SARS-CoV-2 region (yellow box) and low-density SARS-CoV-2 region (white box) in the same plane. Scale bar = 100 µm. All experiments were performed with unique organoid, n = 4 per condition, from the same culturing batch, with images from n = 12 cortical regions with two IPSC lines, and Student's t test was performed (****, P < 0.0001). Experiments were performed twice independently for reproducibility.

Figure 3 .

 3 Figure 3. Single cell RNA-seq of SARS-CoV-2 infected organoids. (A) UMAP projection of cells from single cell sequencing. (B) Monocle trajectory analysis resulted in four distinct states of cells from the organoid. (C) The four major clusters consisted of the following: (1) Neural progenitor, outer radial glia like cells; (2) intermediate progenitor, interneurons; (3) neurons; and (4) cortical neurons. (D) Monocle projection of individual clusters. (E) Heatmap of commonly used genes for identification of cell subtypes in human brain organoids. (F) UMAP projection of organoids depending on infection status. (G) Percentage of infected cells in each cluster. (H) UMAP heatmap of SARS-CoV-2 transcript + cells separated by infection status. (I) Correlation between % change of cell population versus % of cells infected in each cluster. Single cell data were produced by two separate organoids (see F) to ensure reproducibility.

Figure 4 .Figure 5 .

 45 Figure 4. Frequency of each cell cluster from single-cell RNA-seq. Graphs display the percentage frequency of each cluster during a given infection status, 0, 2, 24, and 96 hpi.

  , Fig. 8 B [upper panel], and Fig. 8 D [black arrows]), along with signal in endothelial cells (Fig. 8 B, lower panel).

Amount and quality of data needed: It

  remains one of the main drawbacks to these approaches since the success of DNN critically depends on its training set. However, if classification tasks such as the Image Net challenge required millions of training images to distinguish between hundreds of different classes, segmentation or image reconstruction tasks may require much less but high-quality data.

	Indeed Weigert et al (Weigert et al., 2017) and Wang et al (Wang et al., 2019) show good performances with
	comparatively very little training examples.
	Even though most classical algorithms are outperformed by far by Deep learning (DL) approaches, they still
	present some non-negligible advantages. First, they do not rely on training data and are easily interpretable
	since they are explicit and analytical. Second, they are in a sense more robust where slight tweaks to the input
	data can mislead CNN when not carefully trained.
	Deep learning (DL) still has many challenges and pitfalls to face, far from being the miracle solution to any
	computer vision task. Among the major ones:
	Hallucination problem: Just like the human brain can interpret cloud shapes as animals or faces, Deep
	Neural Networks (DNN) can be tricked by image details when not properly trained. Typically, the presence of
	recurrent undetectable artifacts in the training dataset. GAN for instance is much more sensitive to these
	hallucinations. A strategy to overcome this issue is to add a consistency term to the loss function.
	Overfitting: Overfitting happens when the number of training data is too low. The model memorizes the input
	very precisely and then fails to generalize to unseen data. To overcome this, a typical strategy is to drop out
	regularization to it and to increase the size of the training set. However, this underlies the fact that a too
	homogeneous training dataset can prevent the DNN to perform on slightly different data. Overall a DNN trained
	on a dataset coming from the same specific microscope data might fail to perform very well on data coming
	from another microscope just because it would have slightly different acquisition parameters.

Table 2 :

 2 Summary of the main vasculature extraction pipelines

	features	pipelines Xiang Ji pipeline		Pavone pipeline	VesSAP pipeline (Todorov	CUBIC Cloud (Mano et al.,
			(Ji et al., 2021)				et al., 2020)	2021)
	imaging techniques	serial	two	photon	LSFM imaging	LSFM imaging	LSFM imaging
			tomography
	labeling techniques	lumen perfusion		lumen staining	Vessel walls staining
	imaging resolution	0,3*0,3*1 um3		0,65*0,65*2 um3	1.6 * 1.6 * 3.0
	clearing technique				CLARITY	3DISCO	CUBIC
	segmentation	1 um			0,65*0,65*2 um3	1.6 * 1.6 * 3.0	6,5 um*6,5um*6,5um
	resolution
	sample	deformation	yes
	correction
	scale		whole brain		whole brain	whole brain	Whole mouse
	stitching		Feature-based		Te r a S t i t c h e r	Te r a S t i t c h e r
			registration algorithm
	computational pipeline analytical			Markov random field Deep Learning	(Mano et al., 2021)
			•	feature based
				stitching				•	5 layer CNN	•	machine learning
			•	adaptive
				thresholding			•	registration	voxel
			•	morphological thinning			classification
			•	canny detector	edge			•	ilastix deformation
			•	coherent point drift algorithm			registration
			• •	median filter contrast			•	random forest
				stretching
			•	intensity
				normalization
			•	centerline
				refinement
			•	graph
				correction (self
				loops removal,
				classifier
				based)
			•	registration
	total time		1 brain/week for imaging ?		1 brain/week for imaging	10 days
	Implementation	Matlab2019b		Fiji, Amira 5.3	Matlab2019b
			(Winnubst et al., 2019)	(Winnubst et al., 2019)
	pros		very precise and complete fast		high tissue transparancy,	preserves XFP.
								fast and easy
	cons		long + a bit of an overkill	expensive,	low resolution	Middling transparency.
			maybe, low throughput.	weak	segmentation
			Perfusion based. Tissues	and analysis (Roche
			disruption			et al., 2011)

TABLE d

 d LEAD CONTACT AND MATERIALS AVAILABILITY d EXPERIMENTAL MODEL AND SUBJECT DETAILS

	SUPPLEMENTAL INFORMATION
	Supplemental Information can be found online at https://doi.org/10.1016/j.
	cell.2020.01.028.

B Animals d METHOD DETAILS B Perfusion and tissue processing B Samples staining and iDISCO+ clearing B Blood retention B Ischemic stroke B Light sheet imaging B Computing resources B WobblyStitcher B ClearMap 2.0 toolbox B TubeMap pipeline B Deep vessel filling B Binary smoothing B Skeletonization and graph construction d QUANTIFICATION AND STATISTICAL ANALYSIS B Statistical analysis B Branch point density, branch density, branch length and tortuosity B Shortest paths to arteries and veins B Arterial domains B Distances of vasculature graph features between cortical regions B Stochastic block model analysis d DATA AND CODE AVAILABILITY

  the following we provide a brief summary of the ClearMap 2.0 functionality. ClearMap 2.0 is open source software and fully documented. It is available for download under https://github.com/ChristophKirst/ClearMap2.

	Visualization tools for 3d graphs
	d 3d line plots of 3d graphs
	d 3d mesh plots of 3d vasculature graphs with edge geometries
	Local histogram statistics and 3d rank filter
	Many image processing methods rely on the local histogram or rank order of the voxel intensities in a region around a center voxel
	under consideration. Calculating those histograms for each voxel is computationally expensive and will strongly slow down the pro-
	cessing in images of terabyte size. ClearMap 2.0 thus provides two optimized ways of calculating those local statistics: (i) a rank filter
	toolbox using sliding windows to calculate histograms and (ii) histogram sub-sampling and interpolation routines.
	Alignment methods
	d 3d resampling
	Graph analysis
	d graph preprocessing and cleanup
	d graph branch reduction
	d graph annotation
	d graphs embedded in 3d space
	d 3d graphs with 3d edge geometry
	d 3d edge geometry post-processing, interpolation, and mesh generation
	d morphological operations on edges and vertices
	d sub-graph extraction and spatial slicing
	d network analysis
	On top ClearMap provides a set of visualization tools to visualize large 3d TB datasets and by using overlays or synchronized win-
	dow displays enables direct inspection of the image processing results.
	Visualization tools for 3d images

d 3d alignment to reference atlases (via interface to elastix)

(Klein et al.,

2010) d wobbly stitching (see above) d Allen Brain Atlas annotation modules Image processing methods d clipping and normalization d binary filling d discrete topology based binary smoothing d 3d local gradients and Hessian matrices d 3d tube filter and tubeness measures d 3d rank filter library (> 30 filters) d skeletonization via parallel thinning d 3d tracing d fast calculation of 3d local image statistics d equalization methods d hysteresis and seeded thresholding d 3d adaptive and local image statistics based thresholding d light-sheet artifact removal d fast pseudo deconvolution d deep convolutional neuronal network based 3d image processing d expert processing pipelines for specific applications d fast interactive 2d slice plotting of 3d TB datasets d overlays and/or synchronized window display of multiple datasets. d interactive image processing pipeline generation and parameter exploration GUI d 3d volume rendering d 3d list and line plots e6 Cell 180, 1-16.e1-e13, February 20, 2020

V.3. Evolution of the cerebral vascular density throughout the life of the animal

  One of the goals of my PhD work was to evaluate if changes to neuronal activity levels could impact the structure of the vascular network. This question has not yet received a clear answer from previous studies, with sometimes conflicting results. While it is well demonstrated and accepted that neuronal activity plays a significant role during post-natal development in the construction of the vascular network, it is not yet established if this is the case in a mature, adult brain. The challenge, in answering this question, is that it is difficult to obtain a permanent depression or increase in neuronal activity, as this is a in most cases following homeostatic regulation to prevent the onset of seizures or depressions.

  Downloaded from http://rupress.org/jem/article-pdf/218/3/e20202135/1407835/jem_20202135.pdf by guest on 15 January 2021

Song et al. Journal of Experimental Medicine 5of18 Neuroinvasion of SARS-CoV-2 in humans and mice https://doi.org/10.1084/jem.20202135

  Song et al. Journal of Experimental Medicine 8of18 Neuroinvasion of SARS-CoV-2 in humans and mice https://doi.org/10.1084/jem.20202135 Downloaded from http://rupress.org/jem/article-pdf/218/3/e20202135/1407835/jem_20202135.pdf by guest on 15 January 2021

Please cite this article in press as: Kirst et al., Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature, Cell (2020), https:// doi.org/10.1016/j.cell.2020.01.028

Cell 180, 1-16, February 20, 2020Please cite this article in press as: Kirst et al., Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature, Cell (2020), https://

(legend on next page) 6 Cell 180, 1-16, February 20, 2020 Please cite this article in press as: Kirst et al., Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature, Cell (2020), https:// doi.org/10.1016/j.cell.2020.01.028

Cell 180, 1-16, February 20, 2020Please cite this article in press as: Kirst et al., Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature, Cell (2020), https://

Please cite this article in press as: Kirst et al., Mapping the Fine-Scale Organization and Plasticity of the Brain Vasculature, Cell (2020), https://

Downloaded from http://rupress.org/jem/article-pdf/218/3/e20202135/1407835/jem_20202135.pdf by guest on 15 January 2021

Journal of Experimental Medicine 14 of 18 Neuroinvasion of SARS-CoV-2 in humans and mice https://doi.org/10.1084/jem.20202135 Downloaded from http://rupress.org/jem/article-pdf/218/3/e20202135/1407835/jem_20202135.pdf by guest on 15 January 2021

Journal of Experimental Medicine 15 of 18 Neuroinvasion of SARS-CoV-2 in humans and mice https://doi.org/10.1084/jem.20202135 Downloaded from http://rupress.org/jem/article-pdf/218/3/e20202135/1407835/jem_20202135.pdf by guest on 15 January 2021

ACKNOWLEDGMENTS

The authors are grateful to Carl Modes (MPI Dresden); Marcelo Magnasco and Tobias Neubauer (Rockefeller University); Franca Schmid; Bruno Weber (University of Zurich); Alain Che ´dotal (Vision Institute), Benjamin Charlier, Stanley Durrleman, Laurent Jacob and Jean-Le ´on Thomas (ICM); and Saaid Safieddine and Christine Petit (Institut Pasteur) for insightful discussions. The authors would also like to thank Patricia Gaspar, Silvina Diaz, and Pablo Ariel for their critical reading of the manuscript, as well as Isabelle Bardet, Marie-Claire Camena d'Almeida, and Antoine Renier for their support throughout the project. This work was also made possible by the ICM.Quant imaging core facility (ICM Brain and Spine Institute), in particular with assistance from Basile Gurchenkov and Aymeric Millecamps. This work was supported by the program ''Investissements d'avenir'' ANR-10-IAIHU-06, the ERC starter grant ''NeuroRemod'' (EU Horizon 2020 research program under grant agreement 758817), and the ''Paris Emergences'' program awarded to N.R. C.K. was supported by an independent research fellowship from the Kavli Neural Systems Institute and by an independent physics fellowship from the Center for Physics and Biology at the Rockefeller University. S.S. is supported by the Re ´gion Ile-de-France DIM Math'Innov. A.V.-P. is a Marie Sklodowska Curie Action fellow (Marie Sk1odowska-Curie grant 845685). N.M. is supported by grants from the ANR as part of the second Investissements d'Avenir program LIGHT4DEAF (ANR-15-RHUS-0001) and LabEx LIFESENSES (ANR-10-LABX-65), LHW-Stiftung, and by the Prix Emergence of the Fondation pour l'Audition.

Acknowledgments

We thank the Yale Environmental Health and Safety department for allowing safe working environments with the SARS-CoV-2 virus. We also thank the patient donors, and clinicians who helped with the collection of CSF for neutralization assays. We also thank the members of the Yale Center for Genome Analysis who have helped with all aspects of sequencing. Finally, we thank the members of the Iwasaki laboratory for insightful discussions regarding the project.

This study was supported by National Institutes of Health grants R01AI157488 (A. Iwasaki, S.F. Farhadian), R01NS111242 (A. Iwasaki), T32GM007205 (Medical Scientist Training Program training grant), F30CA239444 (E. Song), 2T32AI007517 (B. Israelow), and K23MH118999 (S.F. Farhadian); the Women's Health Research at Yale University Pilot Project Program (A. Iwasaki and A. Ring); Fast Grant from Emergent Ventures at the Mercatus Center (A. Iwasaki, E. Song, and C.B. Wilen); the Mathers Foundation (A. Ring, C.B. Wilen, and A. Iwasaki); and the Ludwig Family Foundation (A. Iwasaki, A. Ring, and C.B. Wilen). A. Iwasaki is an investigator of the Howard Hughes Medical Institute. Author contributions: E. Song, C. Zhang, K. Bilguvar, and A. Iwasaki planned the project and analyzed data. E. Song, C. Zhang, and A. Iwasaki wrote the manuscript. E. Song, C. Zhang, B. Israelow, A. Lu-Culligan, A.V. Prado, and S. Skriabine performed experiments. P. Lu, O-E. Weizman, F. Liu, Y. Dai, K. Szigeti-Buck, Y. Yasumoto, G. Wang, J. Heltke, E. Ng, J. Wheeler, and M.M. Alfajaro assisted with experiments. E. Levavasseur, S.A.J. Kazmi, K. Zhang, C. Castaldi, B. Fontes, D. Van Dijk, S. Mane, M. Gunel, A. Ring, C.B. Wilen, T.L. Horvath, I. Plu, N.G. Ravindra, S. Haik, J-L. Thomas, A. Huttner, D. Seilhean, A. Louvi, S.F. Farhadian, and K. Bilguvar provided expertise and materials for analysis of data. A.V. Prado, S. Skriabine, and N. Renier performed iDISCO+ imaging and analysis of mice. E. Levavasseur, S.A.J. Kazmi, K. Zhang, I. Plu, S. Haik, J-L. Thomas, A. Huttner, and D. Seilhean provided human samples and analysis of images. A. Iwasaki and K. Bilguvar supervised the project and secured funding.

Disclosures: M. Gunel reported personal fees from AI Therapeutics outside the submitted work; and reported, "AI Therapeutics is currently sponsoring a clinical trial for a therapeutic, which has no relevance for this study, in COVID-19. I am the Chief Scientific Advisor to AI Therapeutics." C.B. Wilen reported personal fees from ZymoResearch outside the submitted work; in addition, C.B. Wilen had a patent for compounds and compositions for treating, ameliorating, and/or preventing SARS-CoV-2 infection and/or complications thereof pending. S. Haik reported a patent to Method for treating prion diseases (PCT/EP 2019/070457) pending. A. Iwasaki reported "other" from RI-GImmune and grants from Spring Discovery during the conduct of the study; in addition, A. Iwasaki had a patent to 14/776,463 pending, a patent for a T cell-based immunotherapy for central nervous system viral infections and tumors pending, and a patent to manipulation of meningeal lymphatic vasculature for brain and CNS tumor therapy pending. No other disclosures were reported.

Submitted: 4 October 2020 Revised: 23 November 2020 Accepted: 10 December 2020

Data availability

Data files from single-cell RNA-seq can be found at Sequence Read Archive accession no. PRJNA682590.

IV. Results

AUTHORS CONTRIBUTIONS

C.K., S.S., A.V.-P., T.T., N.M., M.T.-L., and N.R. designed the study. C.K. wrote TubeMap. S.S. implemented the DNN. A.V.-P. and T.T. performed the experiments. G.G. and P.T. performed the ischemic strokes. C.K., S.S., A.V.-P., P.B., and N.R. analyzed the data. F.V. performed the data annotation. N.R. wrote the initial draft of the manuscript. C.K., S.S., and A.V.-P. contributed sections to the manuscript. All co-authors edited the text.

DECLARATION OF INTERESTS

The authors declare no competing interests.

Supplemental Figures

(legend on next page) Figure S1. Deep Immunolabeling of the Vasculature, Related to Figure 1 (A) Comparison of the tissue and vascular integrities after iDISCO+ processing: coronal 100mm sections of perfused brains (left) or iDISCO+ processed brains (right), immunostained for Neurofilment M and CD31. The iDISCO+ processing of the tissue doesn't visibly affect the vascular or axonal morphology compared to classic immunohistochemistry at this imaging scale. B-E Hemispheres triple immunolabeled with iDISCO+, light sheet scanned and tiled at 4X. A plane is shown close to the midline, as well as an orthogonal reconstructed plane (dotted lines). Insets detail the signal taken at different extremities of the 3D volume, close (i or i') or away (iii or iii') from the imaging and illumination sides. (B) Alternative single labeling of the vessels using rat anti-Podocalyxin only. While this strategy offers the maximal flexibility for counter stains, the signal in the lateral thalamus is very low, but still segmentable (highlighted inset). (C) Alternative double labeling of the vessels against CD31 and Podocalyxin with goat-raised antibodies (D) Alternative triple immunolabeling for the vessels (CD31, Collagen IV and Podocalyxin) further improves the labeling gradient over the dual labeling (CD31, Podocalyxin) and has the best signal to noise ratio, while being less flexible for counter stainings. (E) Validation of the whole-mount passive immunolabeling. Brains whole-mount immunolabeled with the iDISCO+ protocol with primary (Podocalyxin) and secondary rounds of staining, then vibratome sliced and re-stained for CD31 on sections (n = 3). A thalamic section (the most challenging region of the brain in whole-mount stainings), is shown with several insets of different representative regions. Whole-mount and on-section immunolabels overlap perfectly in all regions inspected, surface, deep, as well as in gray and white matter. A rare exception is highlighted: a few capillaries were not detectable in the whole mount staining in the white matter of the cingulate cortex, as shown in E5. Scale bars are 500mm.

(legend on next page) Figure S2. Deep Immunolabeling of Arteries and Veins, Related to Figure 1 (A) Whole brain immunolabeled with iDISCO+ for markers enriched in arteries: sm22 (upper panels) and acta2 (lower panel, alternative host, goat). Lateral projections as well as a 200mm central plane projection are shown. Both markers show comparable staining patterns in the arteries. Expression of both sm22 and acta2 is visible only in a few very large veins. Acta2 expression is also high in the rhinal vein at the cortical surface, while it is sm22 negative. Large veins in general have a patchier distribution and weaker signal than arteries for both markers, which is used as a criteria for later annotations. (B) Overlay of the whole-mount immunolabeling for acta2 (in green, rabbit) with an orthogonal acta2 staining on sections (in red, goat). Both antibodies are raised against the same epitope. A coronal section at the thalamic level is shown. While all arteries are double-labeled, arterioles (low diameter acta2+ vessels) are occasionally only visible from the section staining (arrows, and examples shown in insets). (C) Whole brain iDISCO+ labeled for acta2 and Von Willebrand Factor (vWF). Coronal and sagittal 1mm deep projections are shown. vWF expression is enriched in veins throughout the brain, while also detectable at weaker levels in arteries. The quasi-exclusive presence of veins in the white matter of cerebellar lobules is visible in midline projections. Scale bars are 500mm.

(legend on next page) Figure S3. Maximization of the Imaging Resolution, Related to Figure 1 (A) Evaluation of the signal to noise ratios of capillaries at several depth in a iDISCO+ whole-brain prepared with a triple immunolabeling (CD31, Podocalyxin, Acta2) and scanned at the center of the sample with light sheet microscopy. The images show a constant StN ratio between the surface and the center of the scan at 2500mm. Only at the center of the sample, at 5mm, the scattering starts to be noticeable with a StN ratio of 73. (B) Whole brain stained for Podocalyxin and CD31 imaged with the LaVision LVMI-Fluor 4X/0.3 objective at high light sheet numerical aperture (0.14). Detail of the cross-section of a capillary shown in x,y and x,z reslices, with the corresponding intensity profiles. Depending on its position at the center or on the edge of a tile, the apparent diameter of a capillary will vary between 4 and 6 pixels in the axial direction, but is stable at 3 pixels in the lateral directions. (C) Imaging whole brains with the 12X/0.53 objective generates an apparent diameter of 6 pixels in the lateral directions, closer to the real capillary size, but doesn't improve the axial resolution, which stays unchanged compared with the 4X objective. D-E Tests of the step size to maximize the resolution of the image processing. (D) A tile extracted from a whole brain scan. The maximal intensity projection of the tile, an x,z resliced plane and its binary mask are shown. (E) Detail of the resliced data acquired at different step sizes. Raw resliced data, their binary masks and a mesh projection of the mask are shown. At 2 and 3mm, the resolution is adequate for capillaries, but would fuse the masks of nearby object (red circle). While slightly oversampled for the axial optical resolution, a 1.6mm step-size allows the disambiguation of the masks of some fused vessels. (F) Effect of the light sheet profile on the vascular binary mask. Whole-brain immunolabeling for CD31 and Podocalyxin, iDISCO+ cleared and imaged at a light sheet NA of 0.14. The raw data are presented in the left panel, and the mask on the right side. For each, the original acquisition plane is shown with the (x,z) reslice, highlighting the loss of resolution due to the light sheet heterogeneity along the x-axis. A possible field of view (FOV) crop, between which the binary mask appears isotropic, is suggested (red dotted lines). (G) Effect of the tiling strategy on the number of detected capillaries and their sizes. The tile edges do not introduce a detectable loss of capillary segmented (upper histogram). (H) The diameter distribution along the x-axis shows that most capillaries are measured at a diameter of 5 pixels. The rest of the measured capillary diameters are labeled at 4 or 6 pixels. The detected radii frequency is spatially influenced by the position of the capillary on the tile with a 2 pixels amplitude. As a conclusion, all capillary diameters are measured between 4 and 6 pixels. Scale bars are 50mm (A,B,C) 2 (A) A strategy to generate filled datasets to train the Deep Neural Network: mice are sacrificed unperfused to retain the blood in the vasculature, the skin is removed and the animal immersed in fixative. The brain is dissected after 4 days (middle panel). The blood is clearly visible in the surface veins. The serum is labeled against circulating IgGs, and combined with a podocalyxin co-labeling. (B) Typical results of the whole-blood immunolabeling. Maximum projection through a tile from a whole brain scan. A section plane is indicated in yellow, and details of this plane shown in the center and right panels. The mouse IgG labeling enables the visualization of filled vessels (barely visible in the Podocalyxin labeling, green). The star shows an example of vessel discontinuity in the IgG labeling, corrected from the Podocalyxin labeling. (C) Generation of datasets and model cubes to train a tube-filling neural network: brains are dissected with their vessels still filled with serum through blood retention, iDISCO+ processed and imaged with the same conditions as the regular datasets. The centerlines of the filled vessels from whole brain scans are extracted and used to generate training cubes of artificial hollow vessels, with sizes and densities matching the normal distribution in the brain. (D) Filling of a tube generated from the brain graph (detail from the test cube), and probability map of the added foreground (E) Architecture of the vessel filling DNN. (F) Output of the multipath filter pipeline. A single plane close to the center of the brain is shown. The insets show a caudal region (cerebellum) and rostral (cortex) containing several empty tubes. (G) Output of the tube filling network. The output is based on both the vessel channel and the acta2 channel. Of note, very big veins outside the brain are often missed by the filling, as they do not appear as acta2 positive vessels (an exemple is visible at the ventral surface, close to the circle of Willis). Arteries and smaller veins are efficiently filled by the network. (G) Blocks of raw data annotated by a human annotator (mesh projection of the annotation shown in red). The TubeMap binarization is shown in blue and was performed on the full brain from which the cubes are extracted. Graph are generated from both masks (human and TubeMap), and then the graph branches are compared to the binary mask. Arrows indicate disagreement between the human annotator and TubeMap.

Figure S6

. Automated Annotation of Veins and Artery, Related to Figure 3 (A) Binarization pipeline for the arterial imaging channel (acta2). This binarization pipeline is modeled from the full vasculature binarization (Figure 2A), but adapted to the sparse nature of arterial datasets. (B) Distribution of vessel radii for all vessels (black) or Acta2+ vessels (red). The initial seed points for veins and arteries tracing are placed on large vessels, which are rare in the dataset, and the vessel label is then traced down to smaller radii. (C) Joint distribution of vessel radii and Acta2 signal levels after tracing. Traced veins and arteries are indicated in blue and red respectively (D) Smoothing of the vessel labels can help recover (legend continued on next page) interrupted initial annotations from the Acta2 binary mask (E) Example of the labeling pipeline on a cube of data extracted from a whole brain scan. Raw data are color-coded by depth for clarity. The unlabeled graph is generated from the binary vessel mask. Arterial edges on the full graph are labeled based on their colocalization with the arterial binary mask. Gaps in the initial labeling are then removed by a smoothing operation along the graph edges. Vein label seeds are then added to the large acta2-vessels (initial venous annotation). Smaller arterioles and venules are finally labeled by tracing down the initial annotations to smaller radii.

implementing 3D convolutional DNNs in the laboratory early while trying to solve this problem by extracting the centerline directly via machine learning, instead of topological thinning. I used the previously published "Vessap" network as a starting point, but the results of this network on crossing points were disappointing.

Further work on the design of a robust neural network to detect the centerline of large crossing vessels would be necessary.

Improving the arterial and venous annotation: the last point of necessary improvements to our method is the detection of veins and arteries in the dataset. This is currently done by a complex multi-step pipeline that assigns an arterial identity based on Acta2 (or Sm22) labeling in the raw data, and then labels the veins "negatively" based on the presence of large vessels that do not have strong Acta2 signal. The issue with this strategy is that some veins do express Acta2 at lower levels, and that veins and arteries are not correctly segmented at the brain surface. Therefore, it is complicated to trace an arterial annotation along the graph without incorrectly "jumping" onto a vein. In the initial developments, many vessels had a dual vein-artery annotation.

To i m p r o v e t h i s , w e t e s t e d e x t e n s i v e l y t h e u s e o f a 3 rd positive marker for veins, vWF. While it wasn't possible at the beginning of my PhD to scan with high quality the 3 channels together (CD31, Acta2 and vWF), but this is now possible due to the improvements to the microscope chromatic alignment lenses and scan speeds.

Combining a triple labeling, (a positive marker for arteries and one for veins) with reliable tracing of the annotation along the surface of the brain, may solve the challenge of the arterio-venous labeling. Ideally, it should be possible to seed the annotation of veins and arteries from the brain surface, and then trace the annotation along penetrating vessels until reaching a certain size.

An alternative method could be also to use machine learning on the graph to recognize veins and arteries based on the shape of their vascular trees. Indeed, arteries tend to have branches that make oblique angles from the main tree, while veins have their branches almost perpendicular to the main tree. However, these distinctions are easy to make on the very large veins and arteries but may not apply on smaller vessels.

Validating such approach may also prove difficult.

In conclusion to this part, while the methodology put forward in this work produces high quality reconstructions of the whole brain vasculature, maybe improvements are still possible, and would require much more work soon.

V.2.2 The use of a growth model to understand vascular remodeling

Another strategy to better understand the organization of the vascular network in our data is to design models that simulated the growth of the network, and then use these models to test hypothesis on what factors are driving the local heterogeneities measured, such as variations in vessel densities or orientations.

Most of the models of vascular growth assume the presence of a gradient of signaling molecules, driving the construction of vascular loops in a specific region. One of the most considered hypotheses is that the secretion of these factors is driven by hypoxia for instance. If hypoxia is localized in the tissue, this could constrain the construction of the network. Having the possibility to test in silico different configuration of vessel growth factor gradients would allow us in principle to generate hypothesis on what cell types could be responsible for the plasticity of vascular development.

To t e s t t h i s h y p o t h e s i s , m y g o a l w a s t o c h e c k i f I c o u l d r e p l i c a t e t h e o r g a n i z a t i o n o f t h e v a s c u l a r n e t w o r k i n t h e

sensory cortex, which as a peculiar stratification of vascular features that matches the known layered organization of the neurons. I developed the following model: I randomly seeded "cells" so that the cell density follows a spatial gaussian distribution. Each cell generates a 3D gaussian "VEGF" (ie pro-growth) gradient mimicking its potential requirement for more oxygen. I then generated the corresponding vascular network to evaluate if its topology and physical properties match the ones of a real vascular graph.

In (Lorthois et al 2017), vasculature is modeled as a 3D Voronoi graph built around random seeds evenly distributed in space. This work shows that such a network displays biologically relevant features. I first generated the arterial/venous trees. Seeds are placed at the surface of the space, and then iteratively grown along the Voronoi edges built around the cell's distribution generated earlier. The growth of large vessels along the Voronoi edges follows a simple rule: following the "hypoxic" (or pro-growth) gradient. I then close the network and start the capillary growth process. This is based on simple rules: each vessel is formed by following the Voronoi edges and the gradient. To reproduce the degree 3 distribution of the natural vasculature, each vertex has a probability to fork and become a branch point if its degree is equal to 2. Each vessel tip grows simultaneously until they reach another vessel. The process stops when there are no remaining degree 2 vertices. I then simulated blood flow through the network using the algorithm developed by (Schmidt et al 2018).

Pressure, flow, and velocities are then measured and compared with their counterpart measured from our own data on actual vascular graphs. Results of these comparisons show that the distributions are coherent (Fig 20). Stochastic Block Model partitioning of the model shows that its modularity and number of modules lie within the range observed in real vascular networks. By construction, the graph is 3-regular, and its branch point density matches the cells' distribution. This approach suggests that the gradient-based Voronoi model algorithm may be a biologically relevant vascular model networks and brings improvement to [START_REF] Peyrounette | Multiscale modelling of blood flow in cerebral microcirculation: Details at capillary scale control accuracy at the level of the cortex[END_REF] 's algorithm that didn't force the creation of degree 3 vessels.

ARTICLE

Neuroinvasion of SARS-CoV-2 in human and mouse brain

Eric Song 1 *,CeZh an g 2,3 *, Benjamin Israelow 1,4 , Alice Lu-Culligan 1 , Alba Vieites Prado 5 , Sophie Skriabine 5 ,Pe iw enL u 1 , Orr-El Weizman 1 ,Fe i m eiL iu 1,6 ,Yi l eD ai 1 , Klara Szigeti-Buck 7 , Yuki Yasumoto 7 , Guilin Wang 8 , Christopher Castaldi 8 , Jaime Heltke 8 , Evelyn Ng 8 ,JohnWheeler 8 , Mia Madel Alfajaro 1,9 , Etienne Levavasseur 5 ,BenjaminFontes 10 , Neal G. Ravindra 11,12 , David Van Dijk 11,12 , Shrikant Mane 2,8 , Murat Gunel 2,3,13 , Aaron Ring 1 , Syed A. Jaffar Kazmi 14 , Kai Zhang 14 , Craig B. Wilen 1,9 , Tamas L. Horvath 7 , Isabelle Plu 5,16 , Stephane Haik 5,10,16,17 , Jean-Leon Thomas 5,18 , Angeliki Louvi 3,13 , Shelli F. Farhadian 4,18 , Anita Huttner 19 , Danielle Seilhean 5,16 , Nicolas Renier 5 , Kaya Bilguvar 2,8 , and Akiko Iwasaki 1,15,20  Although COVID-19 is considered to be primarily a respiratory disease, SARS-CoV-2 affects multiple organ systems including the central nervous system (CNS). Yet, there is no consensus on the consequences of CNS infections. Here, we used three independent approaches to probe the capacity of SARS-CoV-2 to infect the brain. First, using human brain organoids, we observed clear evidence of infection with accompanying metabolic changes in infected and neighboring neurons. However, no evidence for type I interferon responses was detected. We demonstrate that neuronal infection can be prevented by blocking ACE2 with antibodies or by administering cerebrospinal fluid from a COVID-19 patient. Second, using mice overexpressing human ACE2, we demonstrate SARS-CoV-2 neuroinvasion in vivo. Finally, in autopsies from patients who died of COVID-19, we detect SARS-CoV-2 in cortical neurons and note pathological features associated with infection with minimal immune cell infiltrates. These results provide evidence for the neuroinvasive capacity of SARS-CoV-2 and an unexpected consequence of direct infection of neurons by SARS-CoV-2.

Single cell-level profiling of human brain organoids

We hypothesized that the cellular heterogeneity of brain organoids may be leading to certain cells being more susceptible to infection and others to death, and that this system would provide an ideal platform to understand the cellular tropism of SARS-CoV-2 in the CNS and elucidate the consequences of SARS-CoV-2 neuroinvasion. To address these questions, we performed single-cell RNA-seq to dissect the cellular states and transcriptional changes occurring after SARS-CoV-2 infection in both infected and noninfected cells within the organoid. We analyzed 60-d-old organoids that were either mock infected or infected with SARS-CoV-2 and collected at 2, 24, and 96 hpi (n = 2 for each time point); 96,205 cells were used to create 31 distinct clusters (Fig. S3 A). Annotations from previous studies [START_REF] Cakir | Engineering of human brain organoids with a functional vascular-like system[END_REF][START_REF] Kanton | Organoid single-cell genomic atlas uncovers human-specific features of brain development[END_REF][START_REF] Velasco | Individual brain organoids reproducibly form cell diversity of the human cerebral cortex[END_REF] allowed us to characterize the identity of each cluster (Fig. S3, A-C). Monocle trajectory further classified these clusters into four major cellular states: neural progenitors/outer radial glia, intermediate progenitor/interneurons, neurons, and cortical neurons (Fig. S3, B-D). To confirm the cell-types we found in the single-cell sequencing, we used additional staining with dorsal cortical markers CTIP2, PAX6, and TBR1 and observed infection primarily overlapping with CTIP2-positive, PAX6-negative cells or CTIP2/TBR1 doublepositive, Pax6-negative cells, which are neuronal cells with a deep-layer fate (layers 5/6). We also observed infection of SATB2positive, CTIP2-negative cells or SATB2/CTIP2 double-positive cells, indicating that upper-layerc o r t i c a ln e u r o n sa r ea l s os u sceptible to infection in the organoid model.

Comparing global UMAP clusters in specimens collected before and after infection, we were able to identify cellular While ACE2 expression levels in the human brain are still being investigated, we showed that ACE2 is expressed at the protein level and is functionally required for SARS-CoV-2 infection in human brain organoids. Further, we detected robust antiviral antibody presence in the CSF of a COVID-19 patient who presented with acute neurological symptoms. This finding suggests that, at least in some patients with COVID-19 and neurological symptoms, there is robust antibody response against the virus within the CSF. In the in vivo setting of the CNS with vasculature and immune cells, neuronal death could have cascading downstream effects in causing and amplifying CNS inflammation.

Although our rodent model does not use endogenous ACE2 expression, it has been previously reported that even mouseadapted SARS-CoV is still neurotropic in wild-type mice, and SARS-CoV-2 is neurotropic in mice with hACE2 expression from the endogenous locus [START_REF] Roberts | A mouseadapted SARS-coronavirus causes disease and mortality in BALB/c mice[END_REF][START_REF] Sun | A Mouse Model of SARS-CoV-2 Infection and Pathogenesis[END_REF]. Using mouse models, we demonstrate for the first time that SARS-CoV-2 neuroinvasion in mice can have significant remodeling of brain vasculature, providing a potential link between the hypoxia and what we see in both the human organoid and the patient brains.

Similar to previous reports of acute hypoxic ischemic damage without microthrombi in postmortem brain of COVID-19 patients [START_REF] Solomon | Neuropathological Features of Covid-19[END_REF], we also found presence of ischemic damage and microinfarcts in postmortem brain samples of COVID-19 patients. In our study, we observed evidence of SARS-CoV-2 infection within the regions of micro-ischemic infarcts, suggesting the possibility of neuroinvasion-associated ischemia and vascular anomalies, consistent with what we observed in mice. However, a limitation of our study is that autopsy samples from only a small number of patients were examined, providing a snapshot of case reports from several patients rather than a generalizable phenomenon. Future studies are needed to examine whether there are other cases of neuroinvasion in the CNS, and the predisposition for such infection. Although we are unable to determine the exact relationship between neuroinvasion and ischemic infarcts, we pose a possible hypothesis from our findings in the patients, mice, and infections of human brain organoids: that SARS-CoV-2 neuroinvasion may cause locally hypoxic regions and disturbance of vasculature, and the disruption of brain vasculature can make vulnerable ischemic infarcts and regions more susceptible to viral invasion (Fig. 8 E). Our findings expand the utility of human brain organoids, beyond modeling fetal brains, and highlight the importance of using a variety of approaches to best model physiology of the human brain.

In future studies, identifying the route of SARS-CoV-2 invasion into the brain, in addition to determining the sequence of infection in different cell types in the CNS, will help validate the temporal relationship between SARS-CoV-2 and ischemic infarcts in patients. It may be through the nasal cavity-to-CNS connection through the cribriform plate, olfactory epithelium and nerve, or viremia, but regardless, the brain should be considered a SARS-CoV-2-susceptible organ system upon respiratory exposure [START_REF] Baig | Potential neuroinvasive pathways of SARS-CoV-2: Deciphering the spectrum of neurological deficit seen in coronavirus disease-2019 (COVID-19)[END_REF][START_REF] Coolen | Early postmortem brain MRI findings in COVID-19 non-survivors[END_REF].

Altogether, our study provides clear demonstration that neurons can become a target of SARS-CoV-2 infection, with devastating consequences of localized ischemia in the brain and cell death, highlighting SARS-CoV-2 neurotropism and guiding rational approaches to treatment of patients with neuronal disorders.

Materials and methods

All procedures were performed in a biosafety level 3 (BSL3)/ animal BSL3 facility (for SARS-CoV-2-related work) with approval from the Yale Environmental Health and Safety committee (#20-19 and #18-16).

planes. (C) Dorsal, ventral, and sagittal projections showing widespread distribution of the virus in the forebrain with patches of high viral density in the cortex (arrow). The virus is not detected in the cerebellum, except for the pial meninges and DCNs. (D) 300-µm-deep projection planes in the cortex showing cortical patches of viral expression (arrowhead), reduced infection of the cells in layer 4, and expression in pyramidal neurons (arrow). (E) ClearMap analysis of the infected cells distribution (n = 3), registered to the Allen Brain Atlas, showing wide distribution of the virus across brain regions, with a few regions with lower densities, among which are dentate gyrus (DG), globus pallidus internal segment (GPi), CA3 hippocampal region, cortical layer 4, and ventromedial hypothalamus (VMH). (F and G) Analysis of the vascular network using ClearMap and iDISCO+ 7 d after intranasal infection by mapping of the vascular network with colabeling of the N protein. Planes at the level of the nose somatosensory cortex are shown. (F) Control uninfected brains. Branch point densities (top panel) peak in controls at layer 4. The density of radially oriented vessels (middle panel) peaks in layers 1, 2, and 3 while decreasing in layers 4, 5, and 6. (G) Brain 7 d after infection. Expression of the N viral protein by neural cells is shown at the level of the nose somatosensory cortex (300-µm projection plane and mapped densities). While branch point densities of vessels still show a peak in layer 4, the normal radial organization of the vessels is not measured in the nose region (arrowhead). Representative render of the vascular graph shows a decrease in vessel orientations in control layers 2 and 3. (H) Schematic of experiment for I and J. Adeno-associated viruses coding for human ACE2 (AAV-hACE2) were injected into the cisterna magna or intratracheally to induce brain-specificorlungspecific expression of hACE2. Brain hACE2-expressing mice were infected with SARS-CoV-2 intraventricularly, and lung hACE2-expressing mice were infected with SARS-CoV-2 intranasally. (I and J) Weight loss curve (I) and survival curve (J) of mice infected with SARS-CoV-2 in the lung (blue) and the brain (red and orange; blue, n = 10; red, n =4;orange,n = 4). Experiments were performed twice for reproducibility. Scale bars = 1 mm (A and C), 200 µm (B), and 500 µm (D and E). CB, cerebellum; DCN, deep cerebellar nuclei; IC, inferior colliculus; SC, superior colliculus; SN, substantia nigra (reticulata or compacta); SS, somatosensory cortex, nose of barrel field.

AAV infection (intratracheal and intracisternal magna injection)

Adeno-associated virus 9 encoding hACE2 (AAV-CMV-hACE2) were purchased from Vector Biolabs.

Intratracheal injection

Animals were anaesthetized using a mixture of ketamine (50 mg/kg) and xylazine (5 mg/kg), injected i.p. The rostral neck was shaved and disinfected. A 5-mm incision was made, the salivary glands were retracted, and the trachea was visualized. Using a 500-µl insulin syringe, a 50-µl bolus injection of 10 11 genome copies (GC) of AAV-CMV-hACE2 was injected into the trachea. The incision was closed with VetBond skin glue. Following intramuscular administration of analgesic (meloxicam and buprenorphine, 1 mg/kg), animals were placed in a heated cage until full recovery.

Intracisternal magna injection

Mice were anesthetized using ketamine and xylazine, and the dorsal neck was shaved and sterilized. A 2-cm incision was made at the base of the skull, and the dorsal neck muscles were separated using forceps. After visualization of the cisterna magna, a Hamilton syringe with a 15°, 33-gauge needle was used to VectaShield mounting medium. Images were acquired using an LSM 880 confocal microscope (Carl Zeiss) and prepared using Fiji.

Electron microscopy

Organoids were fixed using 2.5% glutaraldehyde in 0.1 M phosphate buffer, osmicated in 1% osmium tetroxide, and dehydrated in ethanol. During dehydration, 1% uranyl acetate was added to the 70% ethanol to enhance ultrastructural membrane contrast. After dehydration, the organoids were embedded in Durcupan, and 70-nm sections were cut on a Leica ultramicrotome, collected on Formvar-coated single-slot grids, and imaged on a Tecnai 12 Biotwin electron microscope (FEI).

Single-cell RNA-seq Organoids were collected, and single-cell suspensions were made using a papain dissociation system (Worthington Biochemical Corp.). Single-cell suspensions were loaded onto the Chromium Controller (10x Genomics) for droplet formation. Single-cell RNA-seq libraries were prepared using the Chromium Single Cell 39 Reagent Kit (10x Genomics). Samples were sequenced on the NovaSeq. R v3.4.2 (R Core Team) was used for all statistical analysis. Sequencing results were demultiplexed into Fastq files using Cell Ranger (v3.0.2; 10x Genomics) mkfastq function. Samples were aligned to GRCh38 (10x Genome). The count matrix was generated using the count function with default settings. Matrices were loaded into Seurat v3.1.5 for downstream analysis. Cells with <200 unique molecular identifiers or high mitochondrial content were discarded. Using FindInte-grationAnchors and IntegrateData functions, all libraries were integrated into a single matrix. Principal component values that were statistically significant were identified, and a cutoff point was determined using the inflection point after using the PCElbowPlot function. Clusters were determined using the RunUMAP, FindNeighbors, and FindClusters functions on Seurat. Cells were considered infected if transcripts aligned to viral ORF1ab, Surface glycoprotein (S), ORF3a, Envelope protein (E), Membrane glycoprotein (M), ORF6, ORF7a, ORF8, Nucleocapsid phosphoprotein (N), or ORF10. DEGs from the Find-Markers function were used to perform PANTHER-GO statistical overrepresentation tests for up-regulated and downregulated genes in each condition shown in Fig. 5 B.G e n el i s t s for Fig. 5 C were obtained from gene set enrichment analysis databases (https://gsea-msigdb.org/gsea/index.jsp).

Immunohistochemistry

Paraffin sections were heated for 30 min at 60°C and treated with xylenes, followed by rehydration in decreasing concentrations of ethanol (100, 90, 80, and 70%). Antigen retrieval was performed by boiling in sodium citrate (pH 6.0) for 15 min, and peroxidase activity was blocked with hydrogen peroxide for 10 min. Blocking was performed in 2.5% normal horse serum (Vector Laboratories) and incubated in primary antibody overnight at 4°C. Mouse anti-SARS-CoV-2 spike antibody (clone 1A9; GeneTex GTX632604) was used at a dilution of 1:400. Secondary antibody and detection reagents from the VECTASTAIN Elite ABC-HRP Kit (PK-7200; Vector Laboratories) were used according to manufacturer instructions. Sections were counterstained with Hematoxylin QS (H-3404; Vector Laboratories), dehydrated in increasing concentrations of ethanol, cleared with xylenes, and mounted with VectaMount permanent mounting medium (H-5000; Vector Laboratories).

ACE2 blocking

For ACE2 blocking assays, Organoids were preincubated with anti-ACE2 Antibodies (AF933; R&D Systems) or isotype antibodies at a concentration of 100 µg/ml for 1 h at 4°C. Organoids were then infected with SARS-CoV-2 as described above.

CSF neutralization assay

Excess CSF was obtained from a hospitalized patient with COVID-19 who underwent clinical lumbar puncture and from a healthy control volunteer. CSF was centrifuged to isolate cellfree supernatant, which was used for ELISA and neutralization assays as follows.

ELISA

ELISAs were performed as previously reported [START_REF] Israelow | Mouse model of SARS-CoV-2 reveals inflammatory role of type I interferon signaling[END_REF]. In short, Triton X-100 and RNase A were added to serum samples at final concentrations of 0.5% and 0.5 mg/ml, respectively, and incubated at room temperature for 3 h before use to reduce risk from any potential virus in serum. 96-Well Max-iSorp plates (#442404; Thermo Fisher Scientific) were coated with 50 µl/well of recombinant SARS-CoV-2 S1 protein (100 µg, #S1N-C52H3; ACROBiosystems) at a concentration of 2 µg/ml in PBS and incubated overnight at 4°C. The coating buffer was removed, and plates were incubated for 1 h at room temperature with 200 µl of blocking solution (PBS with 0.1% Tween-20 and 3% m i l kp o w d e r ) .S e r u mw a sd i l u t e d1 : 5 0i nd i l u t i o ns o l u t i o n( P B S with 0.1% Tween-20 and 1% milk powder), and 100 µl of diluted serum was added for 2 h at room temperature. Plates were washed three times with PBS-T (PBS with 0.1% Tween-20), and 50 µl of mouse IgG-specific secondary antibody (1:10,000, #405306; Bi-oLegend) diluted in dilution solution added to each well. After 1 h of incubation at room temperature, plates were washed three times with PBS-T. Samples were developed with 100 µl of TMB Substrate Reagent Set (#555214; BD Biosciences), and the reaction was stopped after 15 min by the addition of 2 N sulfuric acid.

Neutralization assay

Virus for infection was preincubated with 500 µl of healthy or COVID-19 CSF at 37°C for 1 h before infection of organoids. The organoid culture was supplemented with an additional 500 µl of CSF after infection until it was collected for imaging.

Human subjects

Research participants were enrolled at Yale University through Human Investigation Committee Protocols HIC#2000027690 and HIC #1502015318. The Institutional Review Board at Yale approved the protocols, and informed consent was obtained from all participants. Postmortem COVID-19 brain tissues were obtained from COVITIS Biobank (Assistance Publique Hopitaux de Paris, Paris, France). Online supplemental material Fig. S1 shows additional infection of 9-wk organoids with SARS-CoV-2 in addition to infection of hPSC-derived NPCs. Fig. S2 describes in more detail the hACE2 dependence of the virus in infecting neural cells. Fig. S3, Fig. S4,andFig. S5 expand on the findings in human postmortem samples and show more examples of various regions of infected and noninfected brains. Video 1 shows a whole-brain view of SARS-CoV-2 infection.

Supplemental material