**1. Introduction**

The increase of performances of the medical imaging technique in the last decades has allowed non-invasive information of geometries and associated morphologies of large cerebral arteries. In biomedical engineering, this information can be further used for generating 1D to 3D computational models to shed light on cerebrovascular hemodynamics [1,2]. Unfortunately, this process becomes more and more complicated once the vascular scale is

**Citation:** Mañosas, S.; Sanz, A.; Ederra, C.; Urbiola, A.; Rojas-de-Miguel, E.; Ostiz, A.; Cortés-Domínguez, I.; Ramírez, N.; Ortíz-de-Solórzano, C.; Villanueva, A.; et al. An Image-Based Framework for the Analysis of the Murine Microvasculature: From Tissue Clarification to Computational Hemodynamics. *Mathematics* **2022**, *10*, 4593. https://doi.org/10.3390/ math10234593

Academic Editor: Fernando Simoes

Received: 4 November 2022 Accepted: 28 November 2022 Published: 4 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

becoming smaller, reaching the blood–brain barrier capillaries (BBB), whose number is huge (more than 10 billion) [3].

In recent years, extensive research has been oriented toward the microcirculatory flow proposing complex mathematical models based on cerebrovascular images [3–8]. However, image segmentation has inherently several challenges. First of all, patient-specific human images are difficult to be obtained. The acquisition of the images with the necessary resolution in vivo is not feasible at the micro-scale. Furthermore, the use of cadavers for obtaining useful images also affects as the capillaries tend to collapse once the blood flow and pressure reduce after death [9]. For this reason, murine and rodent images have been utilized as baseline geometry for studying BBB microcirculation [10–16]. Other works have presented imaging-driven modeling for hemodynamics in zebrafish microvasculature and mammalian hearts [17–20]. Models oriented to the vascular topology and transport efficiency have been presented by Katifori and coworkers [21,22].

Nowadays, a wide range of imaging techniques are available. In the literature, corrosion casting [23], confocal microscopy [4], computerized tomography angiography and quantitative magnetic resonance angiography [2], two-photon imaging [8,24,25] and synchrotron radiation-based X-ray tomographic microscopy [26–28] have been mostly adopted depending on the specific necessities of the researchers. A combination of some of these methods can thus be of advantage for limiting the weakness of each method and it is applied in the reconstruction protocols of the microvasculature. In this sense, a recent study by Waelchli et al. [29] provides a detailed visualization and quantification of the 3D brain vasculature using resin-based vascular corrosion casting, scanning electron microscopy, synchrotron radiation and desktop microcomputed tomography imaging. These imaging modalities can provide a large field of view of the vascular network but at the same time low resolution. On the contrary, high resolutions are associated with a smaller field [9,30]. Micro-computerized tomography (micro-CT) is a powerful tool for visualizing large vessels but, as aforementioned, it is not capable of imaging properly the microvasculature due to the lack of resolution [31]. For this reason, its use for the microvasculature needs to be modified using additional techniques. With the aim of improving micro-CT performances, Hlushchuk et al. [32] for instance presented an innovative high-resolution micro-CT imaging of animal brain vasculature. Ghanavati et al. [33] proposed a surgical protocol for improving the surgical perfusion of cerebral blood vessels throughout the murine brain and thus obtaining more consistent cerebrovascular images by X-ray micro-CT.

An additional issue is that the cerebral tissue is opaque so that conventional light microscopy is inefficient due to the light scattering provoked by lipids [34,35]. To solve this problem and allow microscopic light to penetrate the brain tissue, several optical techniques have recently been developed. All of them are based on clearing the tissue using chemical procedures. Dodt and coworkers [36] used a mixture of benzyl alcohol and benzyl benzoate to match the refractive index of fixed tissue. However, this protocol only permitted a partial tissue clarification as the clearing solutions led to the rapid loss of fluorescent signals. In a later study, using a so-called 3D Imaging of Solvent Cleared Organs approach (3DISCO), they found that a fast optical clearing can be obtained [37]. Further studies [38,39], using different chemical approaches also achieved rapid tissue clearing encountering similar instability issues.

Tissue-clearing techniques emerged in the last decade to allow high-resolution 3D imaging of biological tissues. Numerous tissue clearing methods are currently available such as DISCO (iDISCO, uDISCO and 3DISCO) [37,40,41], CLARITY [34,42–47], and seeDB [39], among others. Most of these protocols reduce the light scattering provoked by the presence of the lipid and homogenize the RI, obtaining more transparent tissues [48]. Susaki and coworkers [35] developed a whole-brain clearing and imaging method called CUBIC (Clear, Unobstructed Brain Imaging Cocktails and Computational analysis). CUBIC is a comprehensive experimental method involving the immersion of brain samples in chemical cocktails containing aminoalcohols, which enables rapid wholebrain imaging with single-photon excitation microscopy. In parallel, they also improved

their methodology developing the so-called Advanced CUBIC clearing method. This improvement was based on hydration and extended the clearing process to several organs of a mouse, allowing high-resolution 3D imaging [49,50]. Advanced CUBIC was time consuming, had a limited efficiency for clearing organs with high pigment content and adopted the same time for different samples. For these reasons, Res et al. [48] introduced a new ultrasound processing to reduce the clearing time and proposed a new decolorization cocktail to remove pigments. With this optimization method, also called CUBIC-Plus, they enable a considerable shortening of the time acquisition of high-resolution 3D images of the lung. Pinheiro et al. [51] developed an improved clearing protocol, called CUBIC-f, for optimizing fragile samples. Hasegawa et al. [52] introduced CUBIC-kidney for kidney research applications. Based on the CUBIC methodology, Murakami et al. [53] proposed a fluorescent-protein-compatible clearing and homogeneous expansion protocol based on an aqueous chemical solution (CUBIC-X). The expansion of the brain sample allowed the construction of a point-based mouse brain atlas that allows the analysis of numerous samples providing a platform for different organs in the biomedical research, the so-called CUBIC-Atlas [54].

Notwithstanding that imaging techniques are continuously progressing, there are no specific techniques that alone are capable of providing the entire cerebral blood vessels for further reconstruction of comprehensive 3D models [9]. In general, the data obtained after the medical imaging techniques require considerable additional work before these can be treated by computer-aided design programs and computational software. An important pre-processing is for instance necessary for closing all the gaps of the acquired data, simplifying and smoothing the segments of the network that represent the vessels, avoiding noise and generating surfaces that form the limits of the computational domain [3,5,8,9,11,12]. For this reason, in the literature, idealized synthetic computational models based on mathematical algorithms are considered a valuable way to study the cerebral microvasculature. They avoid some of the limitations affecting images-based methods due to the considered microscale [9]. It is about the binary branching trees or networks that mimic the vascular bed morphology. However, the cerebral microvasculature presents for instance loops and anastomoses that cannot be taken into account using simple fractal networks [9]. Hence, simplified binary fractal trees [55–57] have been progressively more and more replaced by complex networks. These models are useful tools for different purposes. In the literature, computer methods have been often based on brain animal images due to the impossibility of invasive experimentation in humans [4,6]. Some studies have been used for interpreting optical measurements acquired in rodents [58,59]. Others were oriented to the analysis of intracellular transport phenomena on length scales not accessible to imaging methods [60,61]. Sherwin et al. [62] introduced a 1D model of a vascular network in space-time variables. Boas et al. [63] presented a symmetric binary vascular network composed by 190 segments to investigate the steady-state or transient response to specific diameter variations of the arteriolar region. Reichold et al. [28] proposed a computational methodology based on anatomical data obtained by synchrotron radiation X-Ray Tomography for simulating rat cerebral blood flow. They presented qualitative results of a fully three-dimensional intra-cortical vasculature structure modeled as a vascular graph. Lorthois and coworkers [3–6,64,65] have provided a large quantitative data focused on the microcirculation of the human cerebral cortex. Recently, they have introduced an analytical model capable of describing the coupling between arteriolar and venular trees, which were modeled using a vascular network approach and the capillary tree, modeled as a continuum porous medium. The research group of Linninger, Hartung and coworkers [2,9,13,14,16,24,66–68] has extensively worked in the microvascular architecture hemodynamics obtaining detailed information on the cerebral microcirculation inside the cerebral cortex by means of vascular networks and numerical algorithms. They analyzed the tissue metabolism coupled to micro-hemodynamics [66], and latterly, they introduced an alternative method to the binary tree for obtaining a more realistic microcirculatory network. This model was generated using Voronoi tessellation first and was later improved

by introducing novel closed networks. It finally includes an arterial and a venous tree with capillary connection synthesized with a single algorithm that allows reducing the computational costs [9].

To the best of our knowledge, there are no studies in the literature that combine tissue clearing, advanced microscopy with images treatment, geometrical reconstruction and numerical simulations. In this study, we aimed to introduce a novel protocol that combines all these techniques for obtaining detailed information about cerebrovascular cortical microcirculation. Concretely, we propose the combination of tissue clearing and advanced microscopy techniques with image treatment, geometrical reconstruction and numerical simulations. With the proposed protocol, we provided a 1D image-based computational model of the cerebral murine microvasculature that allows solving instantaneously the associated hemodynamics. Blood flow features and a quantitative evaluation of the microvascular morphology at the brain cortical territory can be predicted. In particular, different regions and depths of the BBB were considered and investigated with the aim of helping understand its microvascular functionalities and characteristics.

### **2. Materials and Methods**

As described in the previous section, the cerebral tissue is opaque due to the presence of lipids, so conventional light microscopy is inefficient [34,35]. Hence, the lipids need firstly to be removed from this tissue for allowing the light passage without scattering or absorption, matching the refraction index (RI) between the tissue and medium. The tissue clarification is a chemical process of delipidation, decoloring and RI matching. In this work, the mouse brain samples were treated as follows:


Figure 1 shows the step-by-step the process followed. In the next subsections, each step is described.

**Figure 1.** From mouse transcardiac perfusion to staining and cleared tissue.

#### *2.1. Fixation, Sectioning and Tissue Optical Clearing*

Specifically, all animal procedures of this study followed European and Spanish legal regulations and were performed under an ethical protocol approved by the University of Navarra Committee for Ethical Use of Laboratory Animal (076-19). Concretely, *C*57*B*6 mice were euthanized and transcardiacally perfused with PBS (pH = 7.4) and 10 mL of 4% paraformaldehyde in PBS. Mice's brains were dissected, post-fixed overnight in 4% PFA and washed with PBS. Then, 500 μm thick brain slices were sectioned using a vibratome (VT1000S, Leica, Leica Biosystems Technologies, Danaher Corporation, Washington, DC, USA) and kept in PBS solution at 4 ◦C. Sections were cleared following the CUBIC protocol [50]. The method consists of two phases: delipidation with ScaleCUBIC Reagent-1 (urea 25 wt %, Quadrol 25 wt %, Triton X-100 15 wt % and *d*H2O) and refractive index matching with ScaleCUBIC Reagent 2 (urea 25 wt %, sucrose 50 wt %, triethanolamine 10 wt % and *d*H2O). The brain slices were first incubated with ScaleCUBIC Reagent-1 for 4 days at 37 ◦C while gentle shaking, which was followed by PBS washing for 16 h at room temperature. After that, the tissue staining was performed, and the slices were immersed in ScaleCUBIC Reagent 2 until their visualization in the microscope.

### *2.2. Tissue Staining*

Staining of the endothelium in the brain tissue sections was performed between the first and second phases of the clearing process. To this end, delipidated sections immersed in PBS were blocked with BSA 4% and incubated with a 50 μg/mL solution of FITC-Lectin (Sigma-Aldrich, Merck KGaA, Darmstadt, Germany, USA) during 24 h at room temperature. Finally, the sections were incubated with arteriole-specific dye Alexa Fluor 633 hydrazide solution (Thermofisher Scientific, Waltham, MA, USA) (2 μm) for 1 h. Following the arterioles down to the capillaries in the acquired 3D volumes, we determined the areas of connection between the arterial and venous system. This allowed establishing the inlets and outlets of the computational model and hence obtaining later the correct flow direction in the simulations.

### *2.3. Two-Photon Excitation Microscopy*

Image stacks (1 mm × 1 mm × 0.5 mm in size) were collected using a Zeiss LSM 880 (Carl Zeiss, Jena, Germany) equipped with a two-photon femtosecond pulsed laser (MaiTai DeepSee, Spectra-Physics, Milpitas, CA, USA), tuned to a central wavelength of 800 nm, using a 25×/1.8 objective (LD LCI Plan-Apochromat 25×/0.8, Carl Zeiss). Tiles of z-stack scan from 500 μm sections were acquired in the non-descanned mode after spectral separation and emission re-filtering using 500–550 nm and 645–685 nm BP filters for Lectin and Alexa 633 signals, respectively. In Figure 2, a lectin-stained vessels region is shown with a close-up view to a cubic sample.

**Figure 2.** Tissue staining: The tissue was divided in cubic samples of 500 μm: (**a**) endothelial labeling with lectin; (**b**) cubic samples; (**c**) close-up view of the single cubic sample box in red in (**b**).

From the cubic samples represented in Figure 2b (a single region is highlighted in red), five cortical regions were selected and further used for images acquisition, geometrical

reconstruction and computational simulation, as it will be shown in the next sections. The reason why these specific regions were chosen is related to the aim of the study. The cortical regions are areas of the brain located in the cerebral cortex where the BBB is localized and thus are the target of our study because the associated microvasculature is the main filter to pharmacologic drugs. Hence, these cortical regions are of particular interest versus other brain internal areas where there is hardly any capillary network (for example, white matter situated in subcortical regions).

### *2.4. Image Analysis*

The following procedure was applied to each vessel region obtained from the previous steps. All the image treatment was performed in MaTLaB (The MathWorks, Natick, MA, USA), using an appropriate in-house code.

Each volume was first filtered with a 3D Gaussian filter (size: 3 × 3 × 3 voxels, standard deviation: *σ* = 0.5). Even this could lead to a possible loss of information, and the images of the present study are mainly dominated by Poisson noise; thus, we cannot neglect the presence of Gaussian noise that needs to be filtered. Then, a non-local means filter was applied to each slice [69] in the volume to reduce the Poisson noise [70] resulting from the acquisition procedure in the microscope. The vessel-like patterns were enhanced by using morphological filters with linear structuring elements, *Li*, as it was performed in [71], as an adaptation of the top-hat method. In a first stage, an opening is carried out using the previously filtered volume, named as *Sf* , with structuring elements of varying orientation, *Li*. In this way, each *Li* was composed by a length of 51 voxels and a width of 1 voxel. Eight different orientations were defined in the *xy*-plane, i.e., horizontal plane (angular variation was *π*/8, i.e., 22.5◦) to which eight additional angular variations were added in the *z*-axis direction, resulting in 64 possible structuring elements, i.e., *Li* with *i* = 1, . . . , 64.

For each one of the *Li*, a volume was obtained as result of an opening operation. A new volume, named *S*0, was constructed assigning to each voxel the maximum value voxel from the 64 opened volumes previously calculated, as shown in Equation (1):

$$S\_0 = \max\_{i=1,\ldots,64} \left\{ \gamma\_{L\_i} \left( S\_f \right) \right\} \tag{1}$$

where *γ* is the opening operator. This step was finished by means of a geodesic reconstruction using *Sf* as a mask, obtaining the opened volume *Sop*, as shown in Equation (2):

$$S\_{op} = \Gamma\_{S\_f}^{\text{rec}}(S\_0) \tag{2}$$

where Γ is the geodesic reconstruction (opening) operator. The next step was to open *Sf* with the different *Li* and subtract to *Sop*. Then, we added every volume calculated.

$$S\_{\text{sum}} = \sum\_{i=1}^{64} \left( S\_{op} - \gamma\_{L\_i} \left( S\_f \right) \right) \tag{3}$$

Once the vascular structure was enhanced, the binarizing of the volume was performed by using an adaptive threshold [72] of size 71 and a Gaussian statistic (Figure 3a,b). Finally, the biggest connected region was selected and cleaned with a morphological closing (3 × 3 × 3) and a 3D hole-filling approach [73] (Figure 3c).
