2 Introduction

To support healthy brain functioning, the cerebrovascular network undergoes continual adjustments in vessel calibers [1], [2], [3]. Neurovascular coupling refers to the change in blood flow following changes in the level of neuronal activity: under physiological conditions, a generous buffer of nutrients is granted to activated parenchyma via the capillary network [1], [4]. This buffer is maintained through finely tuned regulation of flow through changes in the vessel caliber, mediated via contractile cells in the vessel walls. In the absence of such tuning, pockets of tissue could experience inadequate access to metabolites [5]. Alterations in smooth muscle cells, pericytes, and astrocytes may lead to compromises in vessels’ dilatory capacity and thus deficits in neurovascular coupling [6], [7], [8], [9]. In various brain pathologies, including Alzheimer’s disease, stroke, and trauma, regional blood flow regulation gets impaired through vessel loss and/or dysfunction of the vessels’ dilatory capacity, resulting in regions of ischemia/hypoxia [10], [11]. Previous studies have examined either individual vessels or the tissue level responses, with little attention having been paid to the vascular network, though network dysfunction frequently is associated with accelerated disease progression and long-term symptomatology [3], [12], [13], [14], [15], [16], [17], [18], [19], [20].

While there is copious data on the functioning of individual vessels, interrogation of the microvascular network remains a challenge, in terms of both data acquisition and analysis [6], [7], [21], [22], [23]. To date, studies on the brain vasculature have been done by sparsely imaging individual blood vessels at the cellular scale [6], [7], [21], [24], [25], [26], [27], [28], [29], [30], [31], [32], thereby severely undersampling the microvascular network; or by evaluating the averaged flow over many vessels at the mesoscopic scale, thus failing to discern the flow through individual vessels. A critical gap in the field is the characterization of flow across hundreds of individual vessels, while imaging the network structure that links them together to determine how the vascular response is coordinated across the network. This gap is particularly significant as studies investigating blood flow across several vessels at a time (imaged by varying the line acquisition pattern [6], [24]) have shown highly heterogeneous responses among capillaries. Neuronal function impairments arise wherever local metabolite supply becomes inadequate, notwithstanding the physiological level of flow across the network as a whole, making mapping of vessel changes across the network of particular importance.

To address the limitations of previous work, we developed a novel deep learning (DL)-based pipeline for mapping changes to the geometry of the brain vascular network following neuronal activation, from a time series of volumetric two-photon fluorescence microscopy (2PFM) data. Neuronal activation was elicited by photostimulation of pyramidal neurons expressing Channelrhodopsin-2 (ChR2) [33] in the Thy1-ChR2-YFP mouse model [34]. Our DL pipeline enabled automatic and accurate segmentation, registration and network analysis of large 2PFM datasets across time. We applied our pipeline in a dataset of 17 Thy1-ChR2-YFP mice to map photostimulation-induced changes across the microvascular network - at the level of individual vessels and at the level of vertices spaced every micrometre along the vessels – in relation to the distance to the closest pyramidal neurons expressing the optogenetic actuator and across the cortical depth. Our findings demonstrate the utility of our pipeline for studying in situ microvascular morphology and function to address various neuroscientific hypotheses.

3 Methods

3.1 Cohort

Animals

All experimental procedures in this study followed the ARRIVE 2.0 guidelines [35]. They were approved by the Animal Care Committee of the Sunnybrook Research Institute, which adheres to the policies and guidelines of the Canadian Council on Animal Care and meets all the requirements of the Provincial Statute of Ontario, Animals for Research Act, and the Canadian Federal Health of Animals Act. We used 15 male and 13 female Thy1-ChR2-YFP mice (#007612, line 18, Jackson Laboratory) at 6-12 months of age (283.9 ± 67.0 days), weighing 28.8±7.1g on the day of imaging. The mice were bred in-house and housed under a 12-hr light/dark cycle [34]. In this mouse line, Channelrhodopsin-2 is expressed preponderantly by pyramidal neurons with soma in cortical layer 5 and dendrites projecting to the cortical surface, enabling depolarization of their cellular membranes and action potential generation upon blue light illumination [30], [34]. An attrition schematic is provided in Supplementary Figure 1.

Surgical Preparation

On the day of imaging, mice were induced with 5% isoflurane, transferred to a rectal probe feedback-regulated heating pad (CWE Inc, Ardmore PA) to maintain a temperature of 37 °C, and maintained under 2% isoflurane. A subcutaneous injection of 1 mL of lactated Ringer’s solution was administered at the start of surgery. Throughout the surgical preparation and imaging, we monitored breath rate, heart rate, arterial blood oxygen saturation, pulse distention, and breath distention via a pulse oximeter, with a probe mounted on the thigh (MouseOx, STARR Life Sciences) (Supplementary Table 1). For fine control over respiration, the mice were tracheostomized with an endotracheal tube (20 Ga catheter) and ventilated with a gas mixture of 20-30% O₂ and medical air using a small animal ventilator (SAR 830/P, CWE Inc, Ardmore PA) set to 115-130 breaths per minute at an inspiratory/expiratory ratio of 1:3-4. Their heads were immobilized via ear bars during the placement of the cranial window and imaging. The tail vein was cannulated with a 26 Ga catheter to enable fluorophore and anesthetic delivery. A cranial window was implanted over the forelimb region of the primary somatosensory cortex (AP 0.25mm, ML 2.0 mm); the dura was excised, and 1% agarose was applied between the brain and the glass coverslip. A well was built using dental cement to allow water immersion of the objective. Texas Red dextran (70 kDa MW, Thermo Fisher Scientific Inc, Waltham MA) was diluted in PBS and injected through the tail vein catheter at a concentration of 33 mg/kg for a total injection volume between 0.1 and 0.2 mL. A line containing 0.01 g/ml alpha chloralose was connected to the tail vein catheter. The animal was then positioned underneath the microscope objective, and the structural 2PFM scan acquired, as detailed below. Seventeen mice (nine male and eight female) made it through the surgical preparation (cf. Supplementary Figure 1). Following the structural acquisition, the isoflurane was discontinued, and the continuous infusion of alpha chloralose commenced at 40 mg/kg/hr.

3.2 Imaging

Two-photon Fluorescence Imaging and Optogenetic Stimulation

Mice were imaged on an FVMPE-RS microscope (Olympus, Japan) using a 25x/1.05NA objective. An Insight tunable Ti:Saphire near-infrared laser (SpectraPhysics, USA) was used for 900 nm excitation of Texas Red-labelled vasculature and YFP-labeled ChR2-expressing pyramidal neurons. Two visible light stimulation lasers were used to excite ChR2 at either 458 nm or 552 nm (control). The optical setup is shown in Figure 1. The imaging field-of-view was selected to include at least one penetrating artery and vein yet avoid major pial vessels and thus signal loss in the underlying cortex. Structural scans were acquired under 2% isoflurane using a Galvano scanner, with 2x averaging, a z-step of 0.99 μm, and a nominal lateral resolution of 0.99 μm. Following structural scanning, mice were allowed to rest for 5 minutes to let blood flow equilibrate from brain exposure to room light. The five-minute resting period was observed due to prior reports of red blood cell velocities following optogenetic ChR2 stimulation in the same mouse strain remaining altered for a minute following photostimulation and due to potential pericyte-mediated responses via cytoskeleton reorganization occurring over a minute [6], [30].

Figure 1:

Functional scans were acquired over baseline and post-stimulus periods, flanking a period of blue light illumination. The haemodynamic response to optogenetic stimulation in the Thy1-ChR2-EYFP mouse model was previously observed to last about 45 seconds, and observed to last for minutes in our experimental model (cf. Supplementary Figure 2) [30] These vascular radius time courses were similar to data collected in other optogenetic models of mural cell stimulation where arteries were observed to dilate for approximately 60 and 20 seconds respectively, and capillaries were observed to dilate on the scale of minutes [6], [23]. We thus constrained our volumetric acquisitions to 45 seconds. Each volumetric scan lasted for 42.98 seconds and used the resonant scanner with 5x frame averaging (resulting in a signal-to-noise ratio on the vasculature channel of 6.9±1.6), a nominal lateral resolution of 0.99 μm, a z-step of 2.64 μm, and traversed 250.8 μm of cortical depth. This z-step was chosen because larger z-steps led to discontinuities in the capillary bed segmentation. We imaged the brain in sections from the cortical surface to a depth of 250 μm and from 250 μm to 500 μm of cortical depth. The acquisitions were split into two slabs to maximize the volume coverage without creating discontinuities in the capillary segmentation masks while maintaining a scan duration of under 45 seconds. Stacks were acquired either from the cortical surface down or from the bottom of the stacks toward the cortical surface on different paired acquisitions, so that the post-stimulus delay of image acquisition at different cortical depths was not constant. A 239-μm diameter circular photostimulation region was positioned 250 μm beneath the cortical surface. The photostimulation light was raster scanned over this ROI for 5 seconds between imaging scans, with a pixel dwell time of 4 μs for each pixel of the stimulation ROI. Each pixel in the photostimulation ROI was thus excited every 501.55 milliseconds, i.e. at about 2Hz. We used two laser powers for 458-nm photostimulation: 1.1 mW/mm² and 4.3 mW/mm². The 1.1 mW/mm² is just above the threshold for ChR2 photoactivation and is expected to elicit a small degree of neuronal activity, whereas the 4.3 mW/mm² photoactivation elicits more neuronal activity while still being below ChR2 saturation level [33], [36]. Low powers were used to activate ChR2-expressing neurons to minimize tissue heating [37]. We also performed control experiments via 552-nm photostimulation (where ChR2 excitation falls to 4% of its peak) at 4.3 mW/mm² [38]. Mice rested for approximately 5-7 minutes between scan pairs for the vascular tone to return to its resting state. Photostimulation parameters were presented in randomized order for each mouse.

3.3 Segmentation and Graph Extraction

Ground truth generation

To generate ground truths, we used ilastik’s pixel classification workflow, which relies on random forests, to annotate blood vessels and pyramidal excitatory neurons soma’s in 42 volumes from 25 Thy1-ChR2-EYFP mice for training our DL model (24 images from 16 Thy1-ChR2-EYFP mice from other studies were included in the training, validation, and test cohorts to increase our dataset size). Images were semi-automatically segmented in groups of up to 4 each, with a size of 507×507×250 μm each, due to ilastik’s inability to handle large amounts of annotations over more than a few volumes. The rater labelled targets (neurons vs. blood vessels), corrected mistakes, and classified the pixels where the output was uncertain (as shown using ilastik’s uncertainty guidance feature). During the manual annotation, feature selection was repeatedly optimized using ilastik’s suggest feature function (wrapper method), leading to different features being used for the random forest model for each set of images. Small (under 50 pixels) isolated vascular components were removed with connected component analysis in Python using scikit-image and connected-components-3d. Testing data was withheld until the final model was selected based on the models’ performance evaluation on the validation data set.

Data preprocessing and augmentation

We used the MONAI python package for data preprocessing, augmentation, model training, and prediction [39]. All images and ground truth segmentation masks were up-sampled to an isotropic voxel size using bilinear and nearest neighbour interpolation. Raw 10-bit image intensities were normalized to range from 0 to 1.0. Each volume had eight 128×128×128 pixel patches randomly cropped out for training. Data augmentation, transformations and parameters are listed in Supplementary Table 2. Spatial transformations were selected to expand data variety via cropping, rotations, and mirroring hence exposing the network to images that would be acquired on different positioning of the animal under the microscope. Zooming and deformation transformations were included to expose the network to small changes in the size and morphology of the vasculature. Intensity and Gaussian transformations exposed the network to signal intensity and contrast variations. Dropping pixels was included as the resonant scanner acquisition occasionally yields images with some relatively low signal pixels.

Model Architecture

We trained a state-of-the-art 3D vision transformer (UNETR) model and the U-Net model for baseline comparison, as implemented in PyTorch using the MONAI library [39]. UNETR is a U-Net style architecture with a transformer-based multi-attention head encoder and a CNN decoder [40], [41]. The encoder takes an input image (or patch from each batch) and breaks it down into a sequence of non-overlapping patches, each of size 16×16×16 pixels, which are weighted differently to account for variation in signal intensities and patterns within an image. The sequence of patches is then passed through a multi-head self-attention and multilayer perceptron encoder to capture self-attention between different pixels of the patch to encode long-range relationships between patches. The encoder is then connected to a CNN decoder with skipped connections to the encoder to map features back onto the original image at multiple spatial scales [40]. The model was trained on 2-channel images (EYFP-expressing neurons and Texas Red-labelled vasculature) as inputs to segment neurons and vasculature. We chose the following parameters for the UNETR architecture: a 12 multi-attention head encoder, 16 convolutional features in the first layer of the encoder, a hidden layer size of 768, a multilayer perceptron size of 3072, and Monte Carlo dropout, where on each run of the model 10% of weights were zeroed [40], [41], [42].

The U-net model used is based on the original 2D U-net proposed by Ronneberger et al. in 2015 and extended into 3D via 3D convolutional operations and features residual blocks, parametric rectifying linear units, and instance normalization as described by Kerfoot et al. [41], [43], [44], [45], [46], [47]. The residual blocks were implemented for better gradient flow during the network training. At the same time, the parametric ReLU activation functions were used to improve the ability of the model to adjust its weights during training. Instance normalization was implemented to reduce contrast differences between images fed into the network to improve model robustness. The model was implemented with five layers featuring 16, 32, 64, 128, and 256 channels and dropout. The U-net model was trained using the same data augmentation employed during the UNETR model training.

Model Hyperparameter Optimization, Training, and Prediction

We used a grid search to determine optimal hyperparameters for the UNETR and UNet models within the following parameter space [43], [48]: loss functions (Dice Loss, Dice + Cross Entropy Loss, Dice + Focal Loss, or Tversky Loss), dropout rates (0.1, 0.2, or 0.3), learning rates (5e-3, 1e-3, 5e-4, 1e-4, or 1e-5), and the number of residual units for the Unet model (2 or 3) [49], [50]. We utilized the Adam optimizer [51], and the best model during hyperparameter optimization was selected based on the validation Dice similarity coefficient (DSC) and trained for a maximum of 2400 epochs. The Dice score was the principal evaluation metric, to maximize the overlap between ground truth and prediction masks. Precision and recall were used as secondary metrics to achieve balanced precision and recall where over-segmented results were produced. For early stopping, the epoch with the best performance of the DSC on the validation dataset was selected. The final model that had the best Dice score during hyper-parameter optimization was trained with a Dice + Cross Entropy Loss function, a dropout rate of 0.1, 2 residual units, a learning rate of 1e-5, and a batch size of 1 image with 8 crops per image. Training and optimization were performed on the Narval cluster of Calcul Quebec and the Digital Research Alliance of Canada, with each node using 498 GB of RAM and 4 Nvidia A100 GPUs, each with 40GB HBM2 VRAM.

Ilastik utilized a random forest model from the Vigra library with the default parameters and 100 trees [52], [53]. We initially utilized all default 3D Colour/Intensity/Texture/Edge features during ilastik feature selection and added 2D Colour/Intensity/Texture/Edge features with a sigma of 20.0 pixels. These 2D features were added as a strong, largely planar artifact was observed towards the surface of the images in the neuron channel. We used the wrapper method for feature selection with a set size penalty of 0.10 to determine the optimal features from the starting set. The model was then trained using live updates on the Narval cluster of Calcul Quebec and the Digital Research Alliance of Canada, with each node using 249 GB of RAM and 2 x AMD Rome 7532 CPUs with 64 cores.

For the ilastik model, a subset of 3 training images from 3 mice was used as ilastik’s random forest was unable to train a model with more data even on compute nodes with large memory and CPU resources. The model could never complete training within the time limits of the resource allocations. This inability to complete model optimization was expected as ilastik’s documentation does not recommend training random forest models with full annotation datasets, and increasing the number of annotations does not necessarily lead to better predictions for the pixel classification workflow [52]. The ilastik random forest model was trained with whole images rather than with patches.

Model Comparison

To compare models, we employed several metrics evaluating the similarity between the ground truth and prediction: the Dice Similarity Coefficient (DSC), Precision, and Recall, as well as surface-based metrics: 95% Hausdorff distance (HD95) and mean surface distance. We specifically focused on mean surface distance and HD95 distance since the centerline extraction used in the downstream analysis was highly sensitive to surface irregularities. The model with the lowest standard deviation on the mean surface distance was to be selected absent statistically significant differences in the average mean surface distance between different models.

Consistency in surface placement was key for extracting good centerlines and graphs from the segmentation masks. The metrics utilized are defined as follows:

The DSC measures the overlap between the prediction and ground truth with a value of 1 corresponding to complete overlap, X represents the ground truth and Y represents the predicted value. Precision measures the rate of correctly returned predicted values, whereas recall accesses the rate of return of targeted values. For the Hausdorff 95% distance (HD95), the distance d is the infimum of the Euclidean distance between points x and y (from the surface of segmentation masks X and Y) to segmentation masks Y and X, respectively. HD95 computes the maximum of the 95^th percentile of the supremums (sup) of these minimum surface distances (d(x,Y) and d(X,y)) between the boundaries of the ground truth and the predicted segmentation mask. The mean surface distance measures the average minimum distance between the boundary of the ground truth and the boundary of the predicted segmentation mask. In Equation (5), n_x and n_y are the sets of boundary points of X and Y, respectively. The infimum of the minimum Euclidean distance, d, from boundary points p from n_x, and q from n_y to sets n_x and n_y are then computed and averaged to compute the Mean Surface Distance. Together, the HD95 and Mean Surface Distance assess model performance at the boundaries of the segmentation masks.

Graph extraction of cerebrovascular networks

Graph extraction was performed in Python except for the centerline calculation, which was performed in MatLab (R2021a). Upsampled image acquisitions (0.99 μm isotropic voxel size) were registered with ANTsPy using the rigid registration method with a total sigma of 2 pixels (for smoothing within the registration function) and mean squared error as the similarity metric as input parameters of the registration function (Figure 2A,B). We selected the first baseline scan from each region of interest that was scanned from the bottom to the top to serve as a reference to which all other images from the same region were aligned. The calculated transforms were used to transform images from the same ROI to the space of the reference image using linear interpolation [54]. We then generated segmentation masks for each aligned image using our trained UNETR model with sliding window inference, and we retained Monte Carlo dropout during prediction to create an ensemble of 20 UNETR models, so as to assess the model uncertainty (Figure 2D) [55].

Figure 2.

To extract the centerlines, we assumed that no background was present within the vessels, as regions of background within vessels disrupt the centerline extraction algorithm. Background labelled pixels surrounded by blood vessel segmentations were thus filled in using connected component analysis. Disconnected vascular components under 50 pixels were assumed to be noise and removed. Each aligned binary segmentation mask was dilated three times with a footprint defined by scikit-image’s disk function with a radius of 1 pixel. Next, we computed the union of all segmentation masks from the same ROI generated from registered images (Figure 2E). This common segmentation mask for all time points was eroded to a centerline via a medial axis transformation (bwskel function in Matlab (R2021a)). Critically, using the union of all time points minimized the sensitivity of centreline identification to red blood cell (RBC) stalls at individual time points. “Hair” segments shorter than 20 μm and terminal on one end were removed to reduce false positive vascular branches: this functionality of the skeletonization algorithm is currently unavailable in any python implementation, but is available in Matlab [56]. The vascular centerlines were next used to construct vascular graphs (with sknw) for each ROI, rendering coordinates of vertices along the vessel’s centreline as edges, and branch points as nodes [57].

Based on the vascular graphs, we computed vascular radius estimates at each vertex at each timepoint and then calculated the vertex’s distance to the nearest YFP-expressing pyramidal neuron (as measured by the distance transformation at the vertex to the nearest labelled neuron, with the radius of the vessel subtracted off). We employed a two-stage approach to estimate the radius for each point on the centerline. First, using the binary segmentation mask, we calculated the distance from every vessel pixel in the mask to the nearest background pixel. These values were averaged for each vessel to estimate its radius. The radius estimate defined the size of the Gaussian kernel that was convolved with the image to smooth the vessel: smaller vessels were thus convolved with narrower kernels.

In the second stage, the registered raw image was deconvolved with the point spread function of the microscope, as measured via FluoSpheres carboxylate-modified Microspheres (Cat# F8803, Thermo Fisher Scientific Inc, Waltham MA), using Richardson-Lucy deconvolution. To low pass filter the centreline in advance of computing the tangent vector at each vertex, the coordinates of the vertices along the centerline were smoothed using a Gaussian with a sigma of 3. Next, the tangent vector to the centerline was specified by calculating the gradient of the centerline path. Given this local tangent to the vessel’s direction of travel at a given vertex, we extracted image intensities in the orthogonal plane from the deconvolved raw registered image. A 2D Gaussian kernel with sigma equal to 80% of the estimated vessel-wise radius was used to low-pass filter the extracted orthogonal plane image and find the local signal intensity maximum searching, in 2D, from the center of the image to the radius of 10 pixels from the center. The orthogonal plane image was sampled every 10 degrees (as finer radial sampling did not improve the estimation cf. Supplementary Figure 3) along radial lines emanating from the local signal intensity maximum closest to the center of the image and 5x bicubic upsampled to extract thirty-six 1D signal intensity profiles at that vertex, as shown in Figure 2H. These line profiles were then convolved with a 1D Gaussian kernel with a sigma of 80% of the estimated radius of that vessel, and the gradient of each profile was calculated as shown in Figure 2I, J. The vessel boundary was then placed at the local minimum of the gradient of each profile. The mean and standard deviation of the boundary distances for the thirty-six 1D line profiles were calculated, and boundary points greater than 2 standard deviations away from the mean were excluded (Figure 2K). This radius estimation procedure was repeated for all vertices of all vessels.

In the last step, we computed the distance from each image voxel to the nearest YFP-expressing pyramidal neuron by computing the distance from every image pixel not belonging to the neuron mask to the closest YFP-expressing neuron’s soma boundary (based on the neuron’s binary segmentation mask). At the vessel vertices, these distances were adjusted by subtracting the vessel’s radius. We thereby captured the distance from the vessel surface to the closest YFP-expressing neuron.

Boundary detection validation

To assess the uncertainty in the radius estimates at each vertex, we simulated changes to the vascular diameter by “resizing” the extracted orthogonal plane image and adding Gaussian noise to it. These steps were undertaken to validate the ability of the boundary detection method to estimate diameter changes and evaluate the robustness of the estimates to increased amounts of noise. For the diameter change estimation, images were resampled from the orthogonal plane by a random factor, uniformly distributed in the range of 0.5x to 2x. Then, the aforementioned boundary detection methods were used to estimate vertex-wise calibre changes. The end goal of the assessment was to see how close the boundary detection method-based radius change matched the prescribed diameter change. In the second task, we added Gaussian noise with a sigma randomly chosen from a uniform distribution ranging from 0 to 500 SU to the orthogonal plane image. The noise was added to the image after deconvolution with the PSF and the extraction of the orthogonal plane but before Gaussian smoothing. We then reported on the percent change in the radius as a result of resizing or adding noise in relation to the baseline radius.

3.4 Vascular Network Analysis

Leveraging the graph representation of the vasculature in our pipeline, we next used graph theory to better understand the networks’ behaviour upon neuronal activation. We looked at morphometrics including vascular segment count density, vascular length density, and mean vessel length for comparison with other work. To demonstrate the benefits of extracting a graphical representation of the vasculature in situ, we looked at graph theory metrics including the assortativity of vascular radius changes and changes to the global efficiency of the capillary (below 5 μm in radius) network following optogenetic stimulation. The assortativity measures the tendency for one vessel to dilate or constrict by a similar amount as its neighbours. The assortativity (q) of radius changes in response to stimulation is defined as the Pearson correlation coefficient of these changes on connected vessels:

where e_xy represents the fraction of vessels (edges) in the network that join together nodes with values x and y (i.e. radius changes with values x and y); a_x and b_y are the percentages of edges connecting nodes with values x and y, and σ_a and σ_b are the standard deviations of a_x and b_y [58]. The efficiency, in turn, captures how easily the graph can be traversed. For vascular graphs, efficiency can be conceptualized as the average of the inverse of the total resistive distance of the shortest paths of all combinations of vascular junctions, with the hydraulic resistivity serving as the distance between them. The resistivity (equation 7) is summed along the shortest paths, and the inverse of this sum is then averaged across all shortest paths to compute the efficiency (equation 8). Finally, photostimulation-induced change in the efficiency, from its baseline level, is reported.

In computing resistivity, we assumed a fluidic viscosity (μ) of 4 cP which is within the physiological range [59]; L was the length of the capillary; and R was the radius of the capillary. Vessels greater than 10 μm in diameter were excluded from the efficiency calculation as we wanted to examine blood flow through the capillary nexus between arteries and veins where blood flow may be reversible [60]. The sum of the resistivities along the least resistive path specified ⍴_ij in the efficiency calculation: this sum quantifies how hard it is for blood to move through the said microvascular bed. N refers to the number of nodes in the network. The efficiency was calculated as the average of the inverse of the least resistive paths between all pairs of nodes. The efficiency increases with increasing radius on the shortest paths through the network.

3.5 Statistical models

To statistically compare deep learning model performance, we used the Wilcoxon signed-rank test as implemented in scipy (1.9.3) [61]. Statistics for vascular radius and network metrics were performed in R (4.3.1) using restricted maximum likelihood mixed effects models as implemented in lme4 (1.134) [62], with post hoc comparisons done using emmeans (1.8.8) [63]. We ran the following linear mixed effects model, separately on dilating and constricting vessels that exhibited radius changes above two standard deviations of the vessel’s baseline radius, to examine the effects of photostimulation power on microvascular radius changes in responding vessels:

Due to the difference in their vessel wall ultrastructure, larger microvessels (above 5 μm in radius) were examined separately from capillaries (radius < 5 μm).

For graph metrics, we included nesting of the field of view within each subject as a random effect so as to account for differences in vascular network architecture within an individual. The linear mixed effects models for the graph metrics used were as follows:

In Tables and text, all values have been quoted as mean ± standard deviation, unless otherwise specified.

4 Results

Application of our computational pipeline resulted in robust segmentation of the vasculature and neurons from 4D in situ 2PFM images and rendering of the microvasculature as a graph. The vessel-wise and vertex-wise calibres were tracked across stimulation conditions and related to the cortical depth and the distance to the closest YFP-expressing neuron, mapping the network-level vascular responses to ChR2 activation and revealing the coordination of the microvascular responses following neuronal activation.

4.1 Segmentation model results and comparisons

We compared an ensemble of UNETR models, an ensemble of U-Net models, and an ilastik random forest model on a test dataset of 9 images (507×507×250 μm each) from 6 mice. Examples of segmentation masks produced by each of the models are shown in Figure 4. When evaluating model performance, we paid close attention to the smoothness of the surface of the segmentation masks due to the sensitivity of the centerline extraction algorithms to irregularities in the surface of the masks: smoother vascular segmentation masks resulted in fewer falsely identified vessel branches. Ilastik tended to over-segment vessels, i.e. the model returned numerous false positives, having a high recall (0.89±0.19) but low precision (0.37±0.33) (Figure 3, Supplementary Table 3). When comparing the UNETR and U-Net models, we focused on the surface-based mean surface distance and HD95 distance metrics. Since we observed no significant differences in these metrics between the two models, we selected the UNETR model as our final model because it produced more consistent segmentations on visual inspection and showed significantly better performance than ilastik on HD95 for both vessel and neuron segmentation (p < 0.05).

Figure 3.

Figure 4.

4.2 Vessel extraction improvements via image registration

Rigid registration across all time points from the same field of view improved the ability to trace vascular paths from in situ 2PFM data. Firstly, registration decreased the mean squared error (MSE) between acquisitions from 1306±747 to 0.008±0.003 Signal Units. The number of images acquired per field of view ranged from a total of 2 to 10 depending on how many repeats were able to be acquired. Registration enabled the computation of the union of segmentation masks from all time points. This increased the number of vessel segments identified in each field of view from 241±174 based on a single time point to 412±281 vessel segments per field of view (507×507×250 μm, n=107 fields of view). Taking the union of segmentation masks of an image stack across all time points substantially decreased the incidence of gaps in capillaries, likely arising due to “transient RBC plugs.” The pipeline’s ability to reconstruct the cortical vascular network was thus significantly improved by registering data obtained at different time points.

4.3 Validation of pipeline sensitivity to geometric changes

To evaluate the ability of our computational pipeline to detect vessel caliber changes of various magnitudes, we simulated a range of vascular caliber changes and injected various levels of Gaussian noise. Across >100,000 simulations, the fit of the estimated radius following rescaling against the simulated radius had an R² value of 0.68. Figure 5B presents a heatmap of the estimated radius post-scaling vs. simulated radius, across different vertices of vessel centerlines, highlighting the ability of our pipeline to estimate vascular radii accurately. The addition of Gaussian noise revealed the robustness of the pipeline: radius estimates remained stable with increasing noise levels, until the addition of noise with a standard deviation of over 200 SU (with the intensity of the images ranging from 0 to 1023 SU).

Figure 5.

4.4 Vascular morphology and heterogeneity within and among vessels

Segmentation coupled with graph extraction enabled a detailed characterization of the microvascular network properties. The morphological properties of extracted networks are listed in Table 1, with the probability densities of the vessel length, baseline vessel radius, mean vessel segment depth, and vessel branch point depth shown in Supplementary Figure 4. On the extracted graphs, the vascular radius and distance to labelled neurons were sampled every 1-1.73 μm, enabling detailed analysis of the relationship between the vessel radius change and the proximity to the YFP-expressing neurons. The radius was tracked across different time points, permitting the analysis of the stimulation-induced change in the vascular caliber. To highlight the ability of the pipeline to detect vessels that significantly change their radius after stimulation, Figure 6A shows the standard deviation of the average radii on each vessel segment during baseline frames for three mice. There was a large difference in this standard deviation across various blood vessels, showcasing the model’s ability to reveal baseline variations within each subject. We examined the average change in the vascular radius of each vessel segment after vs. before photostimulation (Figure 6B), with even finer spatial patterns detected by analyzing the vertex-wise radius changes (Figure 6C). The vascular diameter changes were related to the distance from the vessel’s surface to the closest labelled pyramidal neuron at each vertex of the centerline (Figure 6D). The vertex-wise radii estimation allowed the assessment of the variations in radii changes within individual blood vessels (Figure 7). Notably, capillary radius varied along the vessel length across the baseline frames, by 24 ± 28% of the mean resting radius. Consequently, point measurements in vessel calibers - that are widely reported in the literature - do not permit accurate estimation of the microvessel volume changes. Together, the within- and across-vessel radii estimations illustrate the pipeline’s ability to capture spatial variations in the vascular reactivity and relate it to other morphological features (e.g., the distance to the closest labelled neuron).

Figure 6.

Figure 7.

Table 1:

4.5 Vascular reactivity to optogenetic stimulation

The ability of the pipeline to reveal novel spatial relationships between the vascular network reactivity and labelled neurons was demonstrated by examining the relationship between photostimulation-induced microvascular radii responses and 1) the closest YFP-labelled pyramidal neurons within 80μm, and 2) the cortical depth, at the vertex-wise level and across different photostimulations (Figure 8). Vessels were coarsely segregated into small (average radius < 5 μm) and large (average radius > 5 μm) vessels, as we expected them to respond differently due to their differential wall-associated cell composition [6], [64], [65], [66], [67], [68], [69], [70], [71]. Only vessels longer than 20 μm (i.e. vessels whose radius was computed by averaging over many cross-sectional planes) that significantly responded following optogenetic stimulation were analyzed: a vessel was deemed a responder if its radius changed by more than twice the baseline standard deviation in the vessel’s radius. The morphometric properties of the responders, under different stimulation conditions, are listed in Table 2. The average magnitude of significant microvascular radius changes across all stimulation conditions was 1.04 ± 1.11 μm (62±47%). The variability in the radius change within the vessel was higher in the dilating vessels, 0.77 ± 0.61 μm (66±72%) than in the constricting vessels, 0.69 ± 0.49 μm (46±18%), for 458-nm, photostimulation (p<1e-4) (with no statistically significant changes detected for either 458-nm, photostimulation or 552-nm, photostimulation).

Figure 8.

Table 2:

Excluding vessels that did not change their radius by over twice the baseline standard deviation removed almost all large vessels from this analysis. (Notwithstanding, Supplementary Figure 5 depicts unfiltered large vessel constrictions and dilations.) We ran mixed effects models (at the vessel level) separately on constricting and dilating vessels to investigate the effect of optogenetic stimulation power on the vessel radius changes. Each of the models was run separately on small and large vessels.

4.5.1 Vessels further away from labelled neurons constrict while those closer to the activated neurons dilate

We examined the relationship between vascular radius changes and the distance to the closest labelled pyramidal neuron, as microvascular response is thought to result from neuronal activation-elicited generation of vasoactive molecules that diffuse to the neighbouring vessels.

For the control condition (552nm, photostimulation), 2.9% of small capillaries dilated while 1.0% of small capillaries constricted; in larger vessels, barely any responded (0.4% dilated). For this control condition, there was no significant difference in the distance from constrictors or dilators to the closest pyramidal neuron. For the 458 nm photostimulation, capillary constrictors were on average farther away than were dilators from the labelledlabelled pyramidal neuron: dilations occurred 16.8 ± 13.5 μm away from labelledlabelled neurons while constrictions occurred 22.7 ± 16.3 μm for photostimulation (p=1.5e-3) whereas the photostimulation had dilations occur 16.1 ± 14.3 μm away and 21.9 ± 14.6 μm for constrictors (p<1e-4). There was no significant shift between the distance from vessels to neurons for the and stimulations with 458 nm light. Dilations in capillaries following 458 nm photostimulation were larger than those following the 552 nm control photostimulation: 0.90 ± 0.93 μm dilatations occurred with and 0.90 ± 0.78 μm with at 458 nm; vs. 0.58 ± 0.92 μm with at 552 nm (p<1e-4). For constrictions, 458 nm photostimulations led to -1.39 ± 1.51 μm radius changes with and -1.20 ± 1.13 μm radius changes with (p=4.4e-3); whereas 552 nm photostimulation induced -0.37 ± 0.30 μm radius changes with of power, which was smaller than the 458-nm induced responses (p=0.02).

4.5.2 Vascular radius changes at increasing cortical depths

Vascular responses were next segregated by the cortical depth of the vessel (i.e. the average vessel distance from the cortical surface) [72]. Dilators tended to be located closer to the cortical surface across all stimulation conditions. Constricting vessels were located at an average 58 ± 187 μm deeper than dilators for 458-nm stimulation at and 37 ± 179 μm deeper for 458-nm photostimulation at (with no change in the mean depth of either constricting or dilating vessels with changes in the photostimulation power).

4.5.3 Vascular Network Coordination Following Optogenetic Stimulation

We examined the coordination of changes in the microvascular network as a whole via assortativity of radius changes and network efficiency changes. The vessel responses were observed to be assortative, i.e., capillaries mirrored the responses in their neighbours. The increases in stimulation power were accompanied by increases in the assortativity of capillary responses: increasing stimulation level resulted in heightened coordination between adjacent capillaries.

The efficiency increased only at the strongest blue photostimulation, i.e., only at this level of stimulation did the resistivity along the average of all of the shortest paths between junctions in the vascular network decrease, resulting in attenuated resistance to flow through the network. The distribution of changes to the efficiency were highly skewed (with a coefficient of skewness of -1.06 for green illumination, 2.92 for lower intensity blue photostimulation, and 4.87 for higher intensity blue photostimulation). The median increase in the efficiency induced by the higher intensity blue photostimulation, of 4% (IQR: -8% to 38%), was significantly higher than the median -6% (IQR=-9% to 4%) efficiency change following the control green illumination

5 Discussion

Recent studies have demonstrated temporal propagation and coordination in cerebrovascular responses to neuronal activation, whereby arteries dilated after capillary exposure to increased potassium ion concentration [73]; opposing geometric changes have been reported by some studies in capillaries connected by intercapillary tunnelling nanotubes [24]. However, how the effects of these and other mechanisms influence in situ reactivity of the 3D brain capillary network remains unknown. Here, we developed a pipeline for extracting graphs of brain microvascular networks from in situ 2PFM and examine coordination within and across capillaries. Capillary networks and their geometrical changes were imaged via 2PFM during periods of baseline alternated with photostimulation of ChR2 in pyramidal neurons of transgenic mice, and the microvascular mesh was evaluated every 1-1.73 μm. The vascular morphology was then analyzed vertex- and vessel-wise across the entire network. All vessels exhibited significant heterogeneity in caliber changes along their length. Neuronal activation induced both dilations and constrictions of vessels and the incidence of constrictions increased with increasing cortical depth. As the stimulation power increased, the tendency for vessels to change their radius by an amount similar to their neighbours increased. Only the highest photostimulation intensity elicited an increase in the network efficiency. Our findings reveal an intricate level of coordination among brain microvessels and provide a computational analysis platform for interrogating a host of hypotheses on cerebral microvascular reactivity.

Vascular Segmentation and Network Extraction

Intensity thresholding-based image processing pipelines have been used to examine vascular networks and quantify vascular morphology, but they have not gained widespread use due to difficulties in adapting them to highly heterogeneous levels of noise across samples [68], [74], [75], [76]. Deep learning-based methods for analyzing vascular morphology from 3D microscopy images have become prevalent as they provide robust segmentation results over a wide range of signal-to-noise ratios. Recent work has demonstrated steady improvements of segmentation models’ performance with respect to similarity-based metrics (i.e., Dice scores) [77], [78], [79], [80], [81], [82]; though surface-based metrics may be better predictors of how amenable segmentation outputs will be to subsequent morphological analysis of the microvascular network. Our final UNETR model was selected based on a combination of performance metrics including mean surface distance, Hausdorff 95% distance, and rater evaluation, to maximize the smoothness of the surface of generated segmentation masks and reduce false positive branch points during centerline extraction; thereby leading to higher fidelity rendering of microvascular networks. Graph generation was greatly facilitated by computing the union of the vascular segmentation masks across all time points as it enabled tracing of capillaries that had stalls at the individual time points (since the accompanying loss of the fluorescent label otherwise resulted in graph discontinuities).

We tested the ability of the pipeline to detect changes to simulated changes to images and the pipeline’s sensitivity to perturbations to extracted vascular morphology. Using resizing, we confirmed that the boundary detection algorithm was able to detect prescribed changes (cf. Figure 5). The average radius estimates for vessels were then shown to vary by 0.64±3.44% upon changes in centerline position, demonstrating that the radius estimates have a low sensitivity to small (2-3 μm) perturbations in centerline placement. Visual examples of vessels repeatedly dilating or constricting are shown in Supplementary Figures 6 and 7.

Network morphological properties at baseline were in line with prior work. The vascular length density measured from fixed tissue ranged from 0.44 to 1.10 m/mm³ [74], [83], [84], [85], [86], [87]. Our reported vascular density in the forelimb region of the primary somatosensory cortex was 0.40 ± 0.22 m/mm³, with the low end of the range value expected due to fluorescence absorption by hemoglobin in the large pial vessels leading to signal dropout (or shadowing) in the underlying tissue. Our reported average capillary radius of 2.19 ± 1.66 μm was also in line with other studies, where the mean capillary radius ranged from 1.75 to 2.2 μm as measured with confocal microscopy or 2PFM [7], [74].

In Situ Changes to Vessel Calibers Upon Neuronal Activation

Caliber changes at individual vertices along vessel centerlines exhibited significant heterogeneity. Such heterogeneity is expected due to non-uniformly distributed alpha smooth muscle actin-containing cells along vessel walls, as well as differential activations leading to heterogeneous metabolic demand within the tissue [1], [6], [88], [89], [90], [91]. Many previous studies assumed vessel caliber changes to be uniform, compromising the accuracy of the estimates.

As expected, the control 552 nm stimulation led to minimal changes in vessel calibers. To probe for off-target effects, non-transgenic mice were also tested with the same optical setup and photostimulation, with no changes to vascular diameters observed at any of the photostimulation powers utilized (cf. Supplementary Figure 8). In transgenic mice, we detected an average capillary dilation in significantly responding vessels of 70±83% with low-intensity 458 nm stimulation, and 67±61% with higher intensity 458 nm stimulation. Across photostimulation conditions, the capillary dilations ranged from 2% to 805%. These calibre changes were higher than those previously reported, which varied from 2% to 20% depending on the capillary branch order [6], [7], [23], [92], [93]. Far less data are available on constrictions. In the current work, constrictions averaged 47±20% for lower-intensity blue light stimulation and 47±17% for higher-intensity blue light stimulation, with a constriction range of 5% to 97% of the baseline radius. These are higher than the previously reported constrictions of 20% [6], [23], likely due to our identifying as responding vessels only those whose caliber changed by at least twice their baseline caliber variation. It is also worth noting that vessels’ response directions were consistent on repeated trials. Of the vessels whose radius change exceeded twice the baseline variability across time, 31.7% dilated on some trials while constricting on others; 41.1% dilated on each trial; and 27.2% constricted on each trial. (Note that some trials use 1.1 vs. 4.3 mW/mm2 and some have opposite scanning directions).

458-nm photostimulation resulted in a mix of constrictions and dilations with 44.1% of significantly responding vessels within 10 μm of a labelled pyramidal neuron constricting and 55.1% dilating, while 53.3% of vessels further than 30 μm constricted and 46.7% dilated. The cutoff distances from the closest labelled neuron were based on estimates of cerebral metabolic rate of oxygen consumption that showed a steep gradient in oxygen consumption with distance from arteries, CMRO2 being halved by 30 μm away [94]. The stronger blue light stimulation led to an increased rate of constrictions, double that of the low-powered blue light stimulation. For larger vessels, both 458 nm stimulation powers led to a similar dilation level that diminished with increasing distance from labelled pyramidal neurons. This tendency for vessels close to neurons to dilate and further away ones to constrict would be expected in flow redirection into regions of high level of neuronal activity. Stimulation power dependence in blood flow changes has previously been reported in optogenetic mouse models with diffuse stimulation via LED probes, and following transcranial alternating current stimulation [95], [96]. However, neither of the previously employed methods was able to discern the spatial relationship between the vascular caliber changes, or relate these changes to the distribution of the stimulated neurons.

As the blue light stimulation power increased, the mean depth of both constricting and dilating vessels increased, likely resulting from higher intensity light reaching pyramidal neurons deeper in the tissue [97], [98]. The blue light would be expected to excite a lower number of neurons farther from the cortical surface at lower powers. Our results underscore that the haemodynamic response following targeted neuronal activation is not uniformly distributed across the microvascular network: accurate neurovascular coupling assessment thus requires network-based analysis.

Vascular Network Reactivity

To study the microvascular network response as a whole, we examined the assortativity between capillary radius changes and network efficiency changes following optogenetic stimulation. These two graph theory metrics were selected as they both leverage the knowledge of the vascular network structure. Assortativity sheds light on how the vascular network coordinates its responses, while efficiency provides insight into the extent to which those changes facilitate flow through the network. The assortativity revealed that as the stimulation power increased, the tendency of vessels to match their changes to those of their neighbours increased. Previously characterized assortative mechanisms include endothelial cell cation conduction via Kir2.1 channels to synchronize vascular responses [73], and spatial adjacency of pericytes on in vitro retinal preparation leading to assortative changes in neighbouring capillaries [88]. Disassortative (causing opposite changes) mechanisms of capillary coordination have also previously been observed in situ and may result from intercapillary nanotubules’ signalling causing connected pericytes to undergo opposing changes [24]. While not ruling out the presence of disassortative control mechanisms, our results suggest that assortative mechanisms dominate capillary responses to neuronal activation in the somatosensory cortex.

The network efficiency here can be thought of as paralleling mean transit time, i.e., the time it takes blood to traverse the capillary network from the arteries to the veins. In situ studies of mean transit time have revealed a high heterogeneity of plasma traversal of the capillary bed during stimulation, with stimulation reducing plasma transit times by 11% to 20% from its resting levels [93], [99], and simulations suggesting that capillary network geometry and locations of caliber changes exert a substantial influence on these responses [100]. The efficiency of the vascular network here increased significantly only with the strongest 458 nm stimulation. Small dilatations may thus not increase flow in the cortex. The differences in efficiency are likely due to the patterns of localized dilations and constrictions within the vascular network. Efficiency calculations are sensitive to bottlenecks when traversing meshes and certain locations constricting or dilating can have profound impacts on the shortest paths between nodes and the path’s resistivity. The highest powered 458 nm stimulation increasing efficiency may have resulted from increased assortativity causing dilation in key locations within the microvascular network, leading to a significant reduction in shortest path resistivity.

Pipeline Limitations and Adaptability

The segmentation model was trained only on vascular and neuronal labels, limiting its generalizability to segmenting alternative cells in the current state. However, it can easily be fine-tuned or retrained to label other brain cells (e.g., pericytes, astrocytes, or endothelial cells). Our vascular segmentation model generalized well to C57BL/6 mouse and Fischer rat data, as well as to Thy1-ChR2 light-sheet fluorescence microscopy images gathered on an UltraMicroscope Blaze lightsheet fluorescence microscope (Miltenyi Biotech) (cf. Supplementary Figure 9 and Supplementary Table 3). However, the segmentation model performed poorly when significant bleeding occurred in the cranial window, compromising the vascular contrast. Our imaging protocol was challenged by the need to resolve individual vessel responses and capture the entire network within the span of the microvascular response to stimulation: thus, we could not examine the kinetics of the microvascular response. The temporal evolution of the response in individual vessels, however, has been reported on using line scanning acquisitions to measure red blood cell velocity and flux [6], [9], [23], [30], [32], [93]. Moreover, the present analysis pipeline can readily be applied to 2PFM data obtained with finer temporal (e.g., via a Piezo objective positioner) or spatial resolution, and/or different size fields of view. Additionally, alternative definitions of responding vessels may be useful depending on the end goal of a study (e.g. selecting a threshold for the radius change based on a percentage change from the baseline level: cf. Supplementary Figure 10 for capillary changes above 10% of the baseline radius). Finally, microvascular networks in different brain areas may show distinct spatiotemporal profiles of response to neuronal activation. Future work is required to test the generalizability of present findings across different brain regions.

6 Conclusion

We developed a novel deep learning-based computational pipeline for analysis of a time series of 3D 2PFM images and investigation of spatial patterns in microvascular network reactivity to neuronal activation. The microvascular network was represented as a graph, allowing for the evaluation of network geometry changes over time. We tracked the size of blood vessels throughout the network and related vessel radius changes to the distance from the stimulated neurons and the cortical depth. Neuronal activation induced both dilatations and constrictions of capillaries, and the magnitude of these responses increased with increased photostimulation levels while showing significant heterogeneity within and between vessels. In the analysis presented, vertex-wise measurements were aggregated for vessel-wise analysis resulting in highly robust estimates of vessels’ calibers and allowing ready comparisons to literature. Notwithstanding, the pipeline also affords vertex-wise analysis and thus registration of microvascular reactivity with other local morphological features, at an unprecedented spatial scale. With increasing distance of the vessel from the most proximal activated neuron, dilatation magnitude decreased and the incidence of constrictions increased. At the highest stimulation level investigated, the incidence of vessel constrictions also increased with cortical depth. With increasing activation levels, capillaries displayed diameter changes that were similar to their immediate neighbours, while vascular network efficiency increased only under the strongest stimulation. Our computational analysis pipeline permits probing microvascular network reactivity and sheds light on the heterogeneity and coordination of vessel caliber changes across the microvascular network.The pipeline will be made available to the research community to propel future studies of neurovascular coupling and network reactivity.

Author contributions

M.W.R. contributed to experimental conceptualization and design, surgical preparation, imaging and data acquisition, pipeline development, data analysis, and manuscript preparation. J.R.M. contributed to the experimental design conceptualization and manuscript preparation. A.A. contributed to pipeline development and manuscript preparation. A.D. contributed to surgical preparation of mice and manuscript preparation. S.P. contributed to LSFM imaging. M.K. contributed to rat imaging and surgical preparation. M.H. and J.M, contributed to rat imaging. S.P. contributed to LSFM imaging. M.K. contributed to rat imaging and surgical preparation. M.H. and J.M, contributed to Rat imaging. M.G. and B.S. contributed to funding acquisition, experimental conceptualization and design, data analysis, manuscript preparation, and supervised all aspects of the study.

Acknowledgements

We are grateful to Calcul Québec and the Digital Research Alliance of Canada (alliancecan.ca) for their allocation of compute resources that in part supported this research. The authors wish to thank The Imaging Facility, The Hospital for Sick Children, Toronto, Canada for assistance with light sheet fluorescence microscopy used in collecting data. This work was supported by funding from Canadian Institute of Health Research grants PJT376309, PJT156179 and PJT178059. M.W.R was supported by the Queen Elizabeth II/Sunnybrook and Women’s College Health Sciences Centre Graduate Scholarships in Science and Technology. B.S. and M.G. are supported by the Canada Research Chair program (CRC-2018-00042 and CRC-2021-00374).

Data and code availability statement

Data are available upon request, and the code will be made available at the Goubran lab github (https://github.com/AICONSlab).

Conflicts of interest

The authors do not have any conflicts of interest regarding this manuscript, its preparation, or the research that went into it.

Videos

Video 1. Example of vascular segmentation overlaid on raw 2PFM data: https://flip.com/s/_XBs4yVxisNs