Vessel Measurements

Once the volumes were binarized (Figure 3), the main geometrical features were extracted from each region. As first step, the vessel structures were skeletonized using the corresponding morphological operator by obtaining the medial axis [74,75]. In this manner, we could easily define the voxels corresponding to the vessels. Those voxels with only one neighborhood were considered as endpoints, while those others with two neighbors were treated as vessel points. Lastly, those voxels with more than two neighbors were labeled as branchpoints (usually three neighbors, but can be more). The parametrization of our vasculature structure permits breaking our vessel set into individual segments by eliminating the branchpoints that can be more easily analyzed.

**Figure 3.** Image treatment: (**a**) Original images (**left**) and result of filtering (**right**). (**b**) Zoomed versions extracted from (**a**). Original slice (**up**), filtered slice (**down**). (**c**) Segmented volumen. (**d**) Result after artifacts erasing.

Regarding each segment, knowing the voxel size, it was possible to calculate some shape-related characteristics such as the longitude, the curvature or the tortuosity of the segment. Using the vessel set volume and the medial axis voxels data, we also calculated the radius for every skeleton voxel and then computed the mean radius for each segment. Once all the structures were parameterized, we straightforwardly constructed 1D models of the five vessel regions in which each voxel and segmen<sup>t</sup> are characterized for further analyses. A final representative model is depicted in the Figure 3d).

### *2.5. One-Dimensional (1D) Modeling*

### 2.5.1. Governing Equations

The computational models of the previously created geometries were programmed in MaTLaB, and they were based on the hemodynamic network developed by A. R. Pries and T. W. Secomb [76]. In the literature, it is worldwide known that in microvascular networks, the velocities achieved by plasma and red blood cells (RBCs) are very low (60–1 mm/s) [77], translating into a very low Reynolds number (Re < 0.001 for all the analyzed regions) and leading to a laminar capillary flow. Hence, inertial forces have less influence than viscous forces. This fact, in addition with a low Womersley number (*Wo* < 0.01) indicating that the flow can be considered as no pulsatile, enabled a simplification of the Navier–Stokes equation into the Stokes equations (Equation (4)). As a result, we simplified the flow in the capillary bed as a ratio of the pressure drop in every capillary and its hydraulic resistivity:

$$
\mu \nabla^2 \mathbf{v} + \nabla p = 0 \tag{4}
$$

The model handled in this study, after its processing, became a 3D network built from nodes and cylindrical segments in which all the constitutive equations were solved. This model was composed by a collection of interconnected nodes and segments of the BBB microvascular geometry. In this vascular network, the nodes represented locations where vessels bifurcated or ended, and the segments represented the vessels. Each node was defined by coordinates and each segmen<sup>t</sup> was defined by nodes and diameters. Every segmen<sup>t</sup> was then divided into several intermediate segments, making possible the use of tortuous vessel length instead of simplifying the vessels as straight lines between two

end points (Figure 4). Hence, the total length of the vessel was obtained as the sum of the lengths of consecutive intermediate segments (Equation (5), Figure 4).

**Figure 4.** Schematic representation of the segments subdivisions.

$$L\_{ij}^{total} = \sum\_{i}^{k} L\_{ij,k} \tag{5}$$

where *Lij* was the total segmen<sup>t</sup> length divided by *k* intermediate segments. Initially, the mass flow entering a node was the same as the outflow of this node, fulfilling the continuity equation in every geometry node (Equation (6)).

$$\nabla \cdot \mathbf{v} = \mathbf{0} \Rightarrow \mathbf{\dot{m}}\_{\mathrm{in}} = \sum \mathbf{\dot{m}}\_{\mathrm{out}} \tag{6}$$

On the other hand, the blood flow (Q) in every capillary segmen<sup>t</sup> was calculated as shown in Equation (7):

$$Q\_{ij} = \frac{\Delta p\_{ij}}{R\_{ij}}\tag{7}$$

where Δ*p* represented the pressure drop between the defining nodes of the segmen<sup>t</sup> *ij* and *R* represented the flow resistance of segmen<sup>t</sup> *ij*, given a cylindrical shape, which was calculated using the Hagen–Poiseuille Law (Equation (8)):

$$R\_{ij} = \frac{128\mu\_{ij}L\_{ij}}{\pi D\_{ij}^4} \tag{8}$$

The flow resistance of a segmen<sup>t</sup> *Rij* depends on the diameter of the segmen<sup>t</sup> *Dij* and on its length *Lij*. In this study, the tortuous length of each vessel was taken into account instead considering only straight segments adding the tortuosity evaluated during the images' treatment. The flow resistance depends on the blood viscosity of the segmen<sup>t</sup> *μij*, which varies considerably between segments due to the Fahraeus–Lindqvist effect. As known, the latter is caused by the biphasic nature of the blood and the small dimensions of the capillaries [78]. The effective viscosity *μeff* was calculated using the in vivo empirical description made by A. R. Pries [78]. This set of empirical equations takes into account the effects of the biphasic nature of the blood in the capillary bed and calculates its effective viscosity given a vessel diameter *D*, velocity and hematocrit *HD* as follows:

$$
\mu\_{ij} = \mu\_{ij}^{rel} \cdot \mu\_{plasma} \tag{9}
$$

$$\mu\_{ij}^{rel} = \left[1 + (\mu\_{0.45} - 1) \cdot \frac{(1 - H\_D)^{\mathbb{C}} - 1}{(1 - 0.45)^{\mathbb{C}} - 1} \left(\frac{D}{D - 1.1}\right)^2\right] \cdot \left(\frac{D}{D - 1.1}\right)^2 \tag{10}$$

$$
\mu\_{0.45} = 6 \cdot e^{-0.085D} + 3.2 - 2.44 \cdot e^{-0.06D^{0.045}} \tag{11}
$$

$$\mathcal{C} = 0.8 + e^{-0.075D} \cdot \left( -1 + \frac{1}{1 + 10^{-11} \cdot D^{12}} \right) + \frac{1}{1 + 10^{-11} \cdot D^{12}} \tag{12}$$

The hematocrit distribution in the bifurcations of the geometry was calculated using the phase separation law established by A. R. Pries and T. W. Secomb [76]. This law contains a set of empirical equations that define the hematocrit distribution in a bifurcation knowing

the hematocrit in the mother branch and the flow and diameters of the daughter branches (Equations (13)–(15)):

$$FQ\_E = 0, \quad \text{if} \quad FQ\_B \le X\_0 \tag{13}$$

$$\log \text{fit}\_{FQ\_E} = A + B \log \text{it} \left[ \frac{FQ\_B - X\_0}{1 - 2X\_0} \right], \quad \text{if} \quad X\_0 \le FQ\_B \le 1 - X\_0 \tag{14}$$

$$FQ\_E = 1, \quad \text{if} \quad 1 - X\_0 \le FQ\_B \tag{15}$$

where *FQB* was the fractional blood flow in the daughter branch (ratio of the blood flow of the daughter branch and the mother branch) and *FQE* was the fractional erythrocyte flow in the daughter branch (ratio of the erythrocyte flow of the daughter branch and the mother branch). The relationship between erythrocyte flow, blood flow and the segmen<sup>t</sup> hematocrit was obtained using the following equation:

$$Q\_E = Q\_B \cdot H\_D \tag{16}$$

The parameters *A*, *B* and *X*0 were obtained as follows:

$$A = -13.29 \cdot \left(\frac{\frac{D\_A^2}{D\_B^2} - 1}{\frac{D\_A^2}{D\_B^2} + 1}\right) \cdot \frac{1 - H\_D}{D\_F} \tag{17}$$

$$B = 1 + \frac{6.98(1 - H\_D)}{D\_F} \tag{18}$$

$$X\_0 = \frac{0.964(1 - H\_D)}{D\_F} \tag{19}$$

where *DA* and *DB* were the diameters of the daughter branches and *DF* and *HD* were the diameter and hematocrit of the mother branch.

### 2.5.2. Boundary Conditions

It is widely known that one of the most challenging parts in simulating microvascular networks is to establish the conditions in all the in/outflows that appear in the limits of the computational domain. These conditions are necessary for the solution of the 1D equations that describe the blood flow inside the microvasculature. In particular, flow, pressure and hematocrit conditions must be set and, depending on their values, the calculations predict accurate (or less accurate) physiologically meaningful results. In this work, as experimental measurements were not possible in murine brains, some approximations were taken, and literature data were adopted. Different authors used various solutions to this problem [5,24,76,77]. In this work, the solution presented by Lorthois et al. [5] was chosen, as it was simple and fast to implement and achieved valid predictions, comparing the obtained results. The used set of boundary conditions are taken from [77,79] and are described below:

	- 1. At the arterial inflow, a pressure of 50 mmHg was given. The arterial pressure outflow was set to 40 or 45 mmHg depending on its nearness to the inflow. With that, the risk of a short circuit was eliminated.
	- Case 1: Zero flow condition: Flow is set to zero in all the capillary outflows. In this case, the flow goes from the arterial inlet passing through the whole geometry until it reaches a venular outlet. As reported [5], this condition would underestimate the flow in the geometry as it isolates it from its virtual neighbors.
	- Case 2: Constant pressure condition: A constant capillary boundary pressure was calculated so that the net capillary flow (the sum of the flow in all the inlets and outlets) was zero; thus, everything that enters through the arterioles exits through the venules. In other words, this pressure was adjusted such that the total flow entering the arteriolar network was the same as the total flow entering the venular network. In this way, the net flux to all the boundary capillary segments was zero. As a consequence, the net flux leaving the studied brain region through capillaries to supply neighboring areas was exactly compensated by the net flux arriving from neighboring areas through capillaries. As shown in the literature, this condition forces the flux lines to be perpendicular to the ends of the computational domain, maximizing the exchanges of fluid with the neighboring region. For this reason, this condition overestimates the flow in the geometry as it maximizes the flow exchange between the region itself and its virtual neighbors. [64].

The influence of imposing zero flow or a constant pressure at the capillary outlets was found to be limited. Similar results were found also by [5]. Finally, we chose the second option (constant pressure condition), as the first one (zero flow) tended to isolate the volume of the considered capillary regions.
