The mammalian neocortex is innervated by a dense, regulated network of blood vessels known as the cortical angiome . The cortical angiome exhibits neurovascular coupling, namely a temporary change in cerebral blood flow triggered by neuronal activity through direct and indirect signalling pathways, replenishing the surrounding tissue with oxygen and nutrients and removing excess heat and waste[16, 21, 17, 31]. Impaired neurovascular coupling is associated with a variety of debilitating pathological conditions, such as dementia, hypertension, diabetes and Alzheimer’s disease [6, 17, 21]. While the structural properties of the cortical angiome have been considerably elucidated , our understanding of its functional properties is limited, owing in part to the insufficient spatiotemporal resolution of existing imaging techniques . Simply put, it is still unknown how individual microvessels react to individual neuronal action potentials, with preliminary evidence suggesting that the vascular response is mostly driven by specific subtypes of interneurons[29, 9, 21].
Previous attempts to measure the vascular response to neuronal activity were limited to imaging changes in vascular diameter one plane at a time [22, 24, 27, 29, 9]. Most works have repeatedly exposed the animal to an artificial prolonged sensory stimulation over hundreds of consecutive trials, eliciting a vigorous neuronal activation that gave rise to a measurable vascular response [22, 24, 27, 29]. By repeating these visual , olfactory  or somatosensory [27, 29] sensory stimuli while focusing on differing layers of the cortex, a laminar difference in the average onset time and time-to-peak of the vascular response was revealed . In particular, instances of vasodilation begun earlier in the deepest cortical layers, suggesting that vasodilation propagates upwards along penetrating arteries [27, 29, 30].
With neuronal activity per se being at the primary focus of most neuroimaging labs , several imaging methods have been recently tailored for rapidly tracking neuronal activity across considerable large brain volumes, with cellular or near-cellular resolution. These include light field microscopy , lensless imaging , scanned line angular projection microscopy , and reconstruction of 3D imagery from 2D images using deep neuronal networks . While the spatial resolution of these methods is sufficient to discern neuronal cell bodies with little cross-talk, it is insufficient for tracking minuscule changes in cerebral blood diameter, which cause considerable changes in cerebral blood flow. Conversely, while optical coherence tomography and functional ultrasound imaging allow noninvasive, label-free tracking of vascular dynamics over large fields of views, they are incapable of tracking neuronal activity with single cell resolution . The invention of the ultrasonic variofocal lens allowed axial scanning at kHz rates, enabling rapid continuous volumetric multi-photon imaging [19, 13]. It is now technically possible to track the activity of hundreds of neurons and vasoactivity along neighbouring blood vessels, simultaneously across a large brain volume . To name but a few of the benefits of continuous volumetric imaging over traditional planar imaging:
The propagation of vasodilation and vasoconstriction through cortical vessels can be directly observed during spontaneous brain activity, rather than indirectly deduced from averaging repeated evoked trials at differing cortical depths.
Most neuronal cell bodies surrounding a given blood vessel segment are observable with volumetric imaging, but not in planar imaging. Therefore a greater proportion of the neuronal activity that drives vasoactivity is accounted for.
Instances in which neuronal action potentials at a given cortical layer affect metabolic demand at another cortical layer can be accounted for.
Cerebral blood volume can be directly measured for each vessel segment, rather than derived from its diameter along an arbitrary axis in an error-prone fashion .
However, the size of 4D datasets generated by volumetric imaging precludes manual time-lapse segmentation of the imaged blood vessels. Furthermore, the sparse imagery obtained by rapid volumetric two-photon imaging [14, 13] complicates time-lapse vascular segmentation even for trained annotators. The development of accurate algorithms for automated vascular segmentation in 3D-movies is therefore essential for the analysis of neurovascular interactions.
In this paper we show the potential applications of automated angiomal segmentation for the tracking of changes in cerebral blood volume over time, as well as for the identification of spontaneous vascular dilations propagating along penetrating arteries towards the pial surface. To the best of our knowledge, this is the first time that such capabilities are shown.
From a technical perspective, our work extends a recently proposed deep learning technique for microvessel segmentation  called the ACWE network, which is presented in Section 2. This network is trained in an unsupervised manner and does not require carefully labeled training samples, in contrast to other recent deep learning approaches [26, 12, 8, 7, 28]. It is based on the optimization problem minimized by the morphological Active Contours Without Edges method , which is converted into a deep learning solution. The ACWE network was shown to outperform both classical active contour methods as well as the recent deep learning solutions, and to be robust to domain shift across datasets .
While the ACWE network is able to perform well on the task of extracting a single (time-collapsed) microvascular map from a given 4D volume (the same task that is being handled by other recent contributions [11, 12, 26]), there is no prior work capable of handling 4D datasets. We demonstrate that the state of the art static map obtained by ACWE is insufficient for performing an analysis of the temporal dynamics. We therefore propose a method to obtain a sequence of such maps, that makes use of the temporal dynamics, and which is much more suitable for our analysis. The full temporal treatment is made possible by the novel skeleton layer, which ties the segmentation results of the individual frames together.
2 The ACWE network
The ACWE network  reincarnates the ACWE method  as a deep learning technique. This is done by replacing the iterative energy minimization that occurs in the classical method into a loss, and the morphological operations of this method into morphological layers. The network receives an input , which is a 3D intensity-response input volume, where are the volume dimensions of a single intensity channel and the network outputs a segmentation map, , and thresholding is performed to obtain the final result.
The network’s architecture is illustrated in Fig. 1 and consists of a main Encoder-Decoder branch with skip connections, denoted as and (for segmentation), followed by successive operations of Morphological Pooling Layer (Eq. 4-5) for smoothing.
It is trained in an unsupervised manner, using an auxiliary reconstruction loss, provided by an additional decoder , which is used only during training, and outputs a reconstruction .
The network components are then rewritten as follows:
where the operator , is the segmentation before smoothing, and is the segmentation mask obtained after applying the morphological pooling layers and times (the two layers are defined below). The Encoder architecture is based on ResNet34 , where 2D convolutions are replaced with 3D ones. Each of the two decoders and consists of three upsampling blocks with skip connections.
. The layers perform masked max pooling, and then take the maximum or minimum across all results, according to the desired operation ( or respectively). Formally, the function first applies an element-wise multiplication between the mask and the input, denoted by , and then takes the maximum over all locations, see Fig. 3 for an illustration of the 2D case. The layers are define as:
The active contour loss term, , is derived from the ACWE algorithm. Let be the energy that the ACWE minimizes, defined as:
where is the input volume, and () are the average intensities inside and outside the segmentation mask , i.e, , and , where is a voxel. The loss is averaged over all 3D points, where per point in the 3D volume it is given as:
The loss is high if the exponent is applied to a value that is close to zero. This happens if the term is negative and is close to zero, or if is positive and is close to one. Therefore, the loss minimizes . The other loss terms are briefly given as:
Eq. 8 pushes to be significantly larger than , Eq. 9 is the loss of the reconstruction pathway, where is a smoothing term, Eq. 10 pushes the segmentation to be minimal, and Eq. 11- 12 encourage to have semi-binary values close to 0 or close to 1.
The ACWE network employs a compound loss ( are weights):
Our method extends the time-collapsed segmentation , where we add an additional novel skeletonization layer on top of the network and perform per-frame segmentation of the 4D movie. The resulting skeleton serves as an anchoring structure that ties all temporal results together regardless of the transient vascular changes. We first describe our novel differentiable 3D skeletonization layer and proceed with the algorithm specifics.
Skeleton layer The skeleton layer promotes spatially coherent tree-like structures. The novel iterative layer is fully differentiable with respect to the input image, and it is based on the layers and the extension of Lanturjoul’s formula  by Beucher et al.  .
Let be the segmentation output of our network, the skeleton layer’s output is given by , and it is obtained after iterations, starting from :
where for some input , and the operator denotes the open operator, i.e., erosion followed by dilation. The union marks a point-wise maximum, and the erosion operator is denoted by . In Fig. 5(e) we show a sample of the temporal skeleton output.
Training First, network trains on the time-collapsed data obtained by averaging all time frames and generates a segmentation . The skeleton layer is then used () to produce the anchor skeleton from .
The temporal segmentation network is then the same , retrained on each sparse frame of the 4D image (a 3D volume denoted ), with an additional skeleton loss, which encourage the temporal segmentation to be aligned with the static time-collapsed skeleton . For each temporal segmentation made by , we compute the skeleton and the loss:
Pre/Post processing We initially downsample the 4D movie to 60Hz as the input to our network. The resulting 4D segmented movie is then smoothed by a moving average along the z axis with a width of one tenth of the volume.
All imaging experiments and surgical procedures were approved by the Tel Aviv University ethics committee for animal use and welfare and followed pertinent Institutional Animal Care and Use Committee (IACUC) and local guidelines. Please see  for a full description of surgical procedures. Neuronal activity was monitored in an adult male mouse from the C57BL/6J-Tg(Thy1-GCaMP6f)GP5.5Dkim/J transgenic line. Vascular dynamics were monitored by injecting Fluorescein isothiocyanate (FITC) conjugated with 2MDa dextran. Prior to imaging, the mouse was habituated to the imaging conditions across five consecutive days. To minimize motion artifacts during imaging, its head was restrained to a custom-made holder and its platform was clamped to the imaging stage.
We examine two datasets acquired using rapid volumetric two-photon laser scanning microscopy . The first dataset (Fig. 5(a)) tracks cerebral blood volume in nine penetrating arteries and neuronal activity in 103 adjacent neurons (see Fig. 5(a-6)), within a volume of living mouse brain spanning , imaged over 536 sec. at a rate of 125.87 volumes per second. The 4D movie was parsed into voxels and its first 60 sec. (7552 volumes) were selected for analysis. The second dataset (Fig. 5(b)) was acquired within the same imaging session, in the same mouse, at the same magnification, and was centered on the same field of view, albeit spanning only , imaged over 268 seconds at a rate of 125.87 volumes per second. The 4D movie was parsed into voxels and its first 60 seconds (7552 volumes) were selected for analysis. In both datasets, a binning of 2 over time was used for ease of computation, yielding 3776 binned volumes.
All fluorescence values were normalized by their mean and standard deviation, followed by a range stretching, such that the minimal value is 0 and the maximal is 1. A human expert has identified and annotated twelve penetrating arteries in the 3D volume, and we rely on this annotation in our analysis. The human annotation takes the form of rectangular region in each z-slice. We note that our segmentation results could promote the development of automatic vessel annotation methods. However, in this work, we focus on the much more critical bottleneck of obtaining reliable vessel activity measurements.
In Fig. 5 we present the raw data as well as the output of multiple steps of processing. These include: (1) Time-collapsed original data (2) Raw video, (3) Time-collapsed segmentation, (4) Time-varying segmentation (5) Time-varying skeleton. We have also annotated the twelve penetrating vessels. As can be seen, the raw data is almost non informative, while in our time-varying method, the vessels are most clearly visible.
The evaluation is limited by the inability to measure the dynamic behavior by other means. We therefore ran multiple auxiliary experiments as sanity check. In one experiment, we have acquired the same vessels with twice the magnification, in addition to the two 4D movies. To assess our method, we compared the diameter of each annotated vessel, in the temporal segmentation, with its x1 and x2 magnification. For accurate measurements, the ratio in diameters would be exactly 2. With the skeleton layer, the average ratio is . Performing an identical 4D analysis, but without this layer, the ratio is .
|W/O Skeleton layer||0.629||0.585||0.542||0.490||0.871||0.755||0.490||0.307|
|W/ Skeleton layer||0.743||0.731||0.728||0.721||0.956||0.907||0.721||0.485|
Correlation between depth slices While we do not have ground truth vascular data, we can expect certain properties to hold in such data. One easily tested property is the correlation in the measured vascular activity of the same penetrating artery across adjacent axial slices.
As can be seen in Fig. 7(a), the raw measurements are very noisy and show very little correlation between adjacent axial slices. This is further illustrated in Fig. 8(a), which presents the correlation coefficients matrix between different axial slices, after removing the diagonal. The situation is not much better when multiplying the raw data by the segmentation mask of , since the sparsity of the imaging modality leads to very noisy measurements (Fig. 7(b) and Fig. 8(b)). For instance, at a depth of below pial surface, the brightest voxel within the brightest vessel has a time-averaged brightness of merely 0.045 photons per second, corresponding to less than 0.0004 photons per frame in time. As Fig. 7(c) shows, our method leads to a much more coherent dynamic output, which presents a large amount of correlation between adjacent axial slices. Fig. 8(c) shows that this correlation exists only between adjacent axial slices, which is a result of the gradual propagation of waves of vasodilation and vasoconstriction along the penetrating artery.
We tested our proposed skeleton layer for correlation between depth slices. In Fig. 8(c,d) we show the correlation matrix for our method with, and without, the skeleton layer, respectively. As the number of collected photons decreases with imaging depth, the resulting segmentation tends to be less coherent. This is shown in Fig. 8(d) on the bottom right side of the correlation matrix, where the correlation between neighboring slices decreases. Additionally, in Tab. 1, we show the average Pearson correlation for two cases: (i) average correlation as the maximum depth, , increases, considering 5-neighboring slices, and (ii) the average correlation along all depth slices, considering only neighboring slices. In Tab. 2 we show the full experiment results, showing the superiority of the skeleton layer, as an increase percentage in average Pearson correlation, especially in deeper slices and distant neighbors.
Vasodilation and vasoconstriction The dilation and constriction of the penetrating artery are depicted in Fig. 9 and Fig. 7, where we show results for the 12 annotated vessels. As can be seen in panel (a) of Fig. 7, the output of our method, when drawing multiple z-slices on the same plot, varies very fast in time. Indeed our high volumetric sampling rate, used for tracking fast neuronal activity and rigid brain motion, oversampled the propagation of vasoactivity along the penetrating artery. We therefore applied a temporal low-pass filter with a cut-off frequency of 1 Hz. As expected, the low-pass-filtered vasoactivity traces are tightly correlated across axial depths, as shown in panel (b) of Fig. 7. Importantly, instances of vasodilation exhibit an earlier onset at deeper axial slices, as illustrated in Fig. 9(*-2). These observations of individual vasodilations tracked in several axial depths simultaneously are congruent with earlier sensory-evoked observations that were acquired and averaged one plane at a time [27, 29, 30]. Conversely, instances of vasoconstriction exhibit an earlier onset at shallower axial slices, as illustrated in Fig. 9(*-3). As far as we can ascertain, this result was not observed in the past. While Fig. 9 shows only a handful of samples of vasodilation and vasoconstriction, these results are typical, and many more samples were found.
Using automated time-lapse segmentation of blood microvessels in living mouse brain we were able to track, to the best of our knowledge for the first time, how individual instances of vasodilation and vasoconstriction propagate along a penetrating artery. The observed propagation of vasodilation upwards along the penetrating artery is congruent with earlier sensory-evoked observations that were acquired one plane at a time [27, 29, 30]. These propagating waves of vasoactivity along the penetrating artery are not detected by bounding the vessel with a box (Fig. 5(a)), nor by segmenting it using the time-collapsed volume (Fig. 5(b)).
Our ability to track spontaneous vasoactivity along penetrating arteries and other blood vessels in an ecologically-relevant setting paves the path towards linking it with individual action potentials. By computing the spike-triggered vasoactivity for different neuronal subtypes, we plan to distill a canonical hemodynamic response function (HRF), namely the small-signal impulse response of various classes of vessels to a single neuronal action potential.
This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant ERC CoG 725974). PB acknowledges the ERC (grant 639416) and the Israel Science Foundation (grant 1019/15) for support of different aspects of this project. The authors thank David Kain for conducting the mouse surgery. SG’s contribution is part of a Ph.D. thesis research conducted at TAU.
-  (2012) The brain activity map project and the challenge of functional connectomics. Neuron 74 (6). Cited by: §1.
-  (2018) DiffuserCam: lensless single-exposure 3d imaging. Optica 5 (1), pp. 1–9. Cited by: §1.
-  (1994) Digital skeletons in euclidean and geodesic spaces. Signal Processing. Cited by: §3.
-  (2013) The cortical angiome: an interconnected vascular network with noncolumnar patterns of blood flow. Nature neuroscience. Cited by: §1.
-  (2001) Active contours without edges. IEEE Trans. image processing. Cited by: §1, §2, §2.
-  (2018) The effect of hyperglycemia on neurovascular coupling and cerebrovascular patterning in zebrafish. J Cerebral Blood Flow&Metabolism. Cited by: §1.
-  (2018) Automatic graph-based modeling of brain microvessels captured with two-photon microscopy. j. biomedical and health informatics. Cited by: §1.
-  (2018) Whole-brain vasculature reconstruction at the single capillary level. Scientific reports. Cited by: §1.
-  (2019) An oligarchy of no-producing interneurons controls basal and evoked blood flow in the cortex. bioRxiv, pp. 555151. Cited by: §1, §1.
-  (2014) Determination of vessel cross-sectional area by thresholding in radon space. Journal of Cerebral Blood Flow & Metabolism 34 (7). Cited by: item 5.
Unsupervised microvascular image segmentation using an active contours mimicking neural network. In
IEEE International Conference on Computer Vision (ICCV), Note: Available online as an arXiv preprint 1908.01373. Cited by: §1, §1, §2, Figure 7, §3, §4.
Deep convolutional neural networks for segmenting 3d in vivo multiphoton images of vasculature in alzheimer disease mouse models. arXiv preprint 1801.00880. Cited by: §1, §1.
-  (2019) Improving in vivo multi-photon microscopy using plug and play photon counting. pp. JT4A–10. Cited by: §1, §1.
-  (2018) PySight: plug and play photon counting for fast continuous volumetric intravital microscopy. Optica 5 (9). Cited by: §1, §4, §4.
-  (2016) Deep residual learning for image recognition. CVPR. Cited by: §2.
-  (2018) What is the key mediator of the neurovascular coupling response?. Neuroscience & Biobehavioral Reviews. Cited by: §1.
-  (2019) Neurovascular and cognitive dysfunction in hypertension: epidemiology, pathobiology, and treatment. Circulation Research. Cited by: §1.
-  (2019) Kilohertz frame-rate two-photon tomography. Nature Methods. Cited by: §1.
-  (2015) Continuous volumetric imaging via an optical phase-locked ultrasound lens. Nature methods. Cited by: item 4, §1.
-  (2010) On the use of the geodesic metric in image analysis. Journal of Microscopy. Cited by: §3.
-  (2019) How reliable is cerebral blood flow to map changes in neuronal activity?. Autonomic Neuroscience. Cited by: §1.
-  (2016) Neural correlates of single-vessel haemodynamic responses in vivo. Nature. Cited by: §1.
-  (2016) Compressive light-field microscopy for 3d neural activity recording. Optica. Cited by: §1.
Vascular compartmentalization of functional hyperemia from the synapse to the pia. Neuron. Cited by: §1.
-  (2019) Computational processing of neural recordings from calcium imaging data. Current opinion in neurobiology 55, pp. 22–31. Cited by: item 4.
-  (2016) Deep learning convolutional networks for multiphoton microscopy vasculature segmentation. arXiv preprint arXiv:1606.02382. Cited by: §1, §1.
-  (2010) Cortical depth-specific microvascular dilation underlies laminar differences in blood oxygenation level-dependent functional mri signal. PNAS. Cited by: §1, Figure 9, §4, §5.
Automated analysis of whole brain vasculature using machine learning. bioRxiv, pp. 613257. Cited by: §1.
-  (2016) Cell type specificity of neurovascular coupling in cerebral cortex. Elife. Cited by: §1, §1, Figure 9, §4, §5.
The roadmap for estimation of cell-type-specific neuronal activity from non-invasive measurements. Phil. Trans. of the Royal Society B: Biological Sciences. Cited by: §1, Figure 9, §4, §5.
-  (2017) Understanding the neurovascular unit at multiple scales: advantages and limitations of multi-photon and functional ultrasound imaging. Adv. drug delivery reviews. Cited by: §1, §1.
-  (2019) Three-dimensional propagation and time-reversal of fluorescence images. arXiv preprint arXiv:1901.11252. Cited by: §1.