Brain Capillary Networks Across Species: A few Simple Organizational Requirements Are Sufficient to Reproduce Both Structure and Function

Despite the key role of the capillaries in neurovascular function, a thorough characterization of cerebral capillary network properties is currently lacking. Here, we define a range of metrics (geometrical, topological, flow, mass transfer, and robustness) for quantification of structural differences between brain areas, organs, species, or patient populations and, in parallel, digitally generate synthetic networks that replicate the key organizational features of anatomical networks (isotropy, connectedness, space-filling nature, convexity of tissue domains, characteristic size). To reach these objectives, we first construct a database of the defined metrics for healthy capillary networks obtained from imaging of mouse and human brains. Results show that anatomical networks are topologically equivalent between the two species and that geometrical metrics only differ in scaling. Based on these results, we then devise a method which employs constrained Voronoi diagrams to generate 3D model synthetic cerebral capillary networks that are locally randomized but homogeneous at the network-scale. With appropriate choice of scaling, these networks have equivalent properties to the anatomical data, demonstrated by comparison of the defined metrics. The ability to synthetically replicate cerebral capillary networks opens a broad range of applications, ranging from systematic computational studies of structure-function relationships in healthy capillary networks to detailed analysis of pathological structural degeneration, or even to the development of templates for fabrication of 3D biomimetic vascular networks embedded in tissue-engineered constructs.


S1.1 Generation of Voronoi diagram
Each cell of length L C was divided into 16 × 16 × 16 pixels, and one seed was placed at a random pixel. The 3D Voronoi diagram was extracted using the MATLAB function voronoin, which returns the coordinates of vertices and the list of vertices belonging to each polyhedron, from which the network connectivity matrix was constructed.

S1.2 Pruning the network
Very small or narrow polyhedral faces (with area < 75 pixels or with any interior angle < 30°) were merged with the neighboring face sharing the longest edge of the current face ie. this edge was deleted. Similarly, neighbouring faces lying almost in a plane (solid angle < 15°) were merged. If any triangles with area < 75 pixels remained, the longest edge was removed. Final networks were not very sensitive to the choice of face area or angle criteria. Vertices were then removed or merged in two stages: 1. Boundary vertices. If two boundary vertices v b, 1 and v b,2 were located within some minimum distance D min of each other, v b,1 was arbitrarily deleted, but only if it was connected to an interior (i.e. non-boundary) vertex v i,1 which had more than three connections. If both boundary vertices were connected to the same interior vertex v i,1 , which itself only had three connections (v b,1 and v b,2 plus one other interior vertex v i,2 ), then v b,1 and v i,1 were deleted, and v b,2 was instead connected to v i,2 . If neither of these cases applied, v b,1 was moved a distance D min along the boundary away from v b,2 in the direction of the vector between the two vertices.
2. Interior vertices. To minimize the computationally heavy task of calculating the distance between all vertices, the domain was divided into sub-cubes of size L C × L C × L C (solid lines in 2D in Figure S1a). These cubes were shifted by 0.5L C in each direction yielding an off-set grid, to avoid missing vertex pairs spanning neighboring cubes (dashed lines in Figure S1a). Running through cubes in randomized order, the distance between each pair of vertices within the same cube was calculated. If a pair of vertices (v 1 and v 2 ) lay < D min from each other, as in Figure S1b, v 1 was arbitrarily deleted and its connections were re-attached to v 2 ( Figure S1c), but only if v 2 still had at least three connections. Otherwise, v 1 was moved a distance D min away from v 2 in the direction of the vector between the two vertices.
Final networks were not very sensitive to the choice of D min , and thus its value was fixed at 2.5 pixels (for L C = 75µm, D min = 11.7µm). Next, excess edges between internal vertices were removed in random order according to the following criteria. An interior vertex could not have more than one connecting edge with a boundary vertex, except if it only had one internal connection, in which case it may be connected to two boundary vertices. If an internal vertex was not connected to any boundary vertices, it should have three connections to internal vertices.

S1.3 Splitting multiply-connected vertices
Vertices with more than 3 connections were treated in random order. For each multiply-connected vertex, the shortest connected edge was kept, as well as the edge with smallest solid angle relative to that Figure S1: a) Schematic diagram illustrating the division of the domain into a grid of length L C (solid lines) and an off-set grid shifted by 0.5L C (dashed lines), for efficient detection of vertices (black dots) less than distance D min apart. b) Illustration of two vertices v 1 and v 2 separated by less than D min ; c) vertex v 1 and its connections are merged into vertex v 2 ; d) -g) Sequence of schematic diagrams illustrating the division of a multiply-connected vertex with five connections into three bifurcations. See text for further details.
edge (see e.g. the two edges in bold in Figure S1d). The center of mass of the remaining vertices was calculated (the cross in blue in Figure S1e), and a new vertex was placed on the vector between the initial vertex and the center of mass, at a distance of half the shortest length of the remaining edges, or D min , whichever is smaller (see the blue arrow in Figure S1e). This process was repeated until all interior vertices had only three connections (Figure S1f and g). This step introduced more short vessel lengths, but since there were typically only a small percentage of multiply-connected vertices, this did not have an important effect on overall network properties.

S1.4 Domain size
To avoid any boundary effects associated with the generation of these synthetic networks, for a desired domain size of L x × L y × L z , a larger network was initially generated in a domain of 2(L x + 2L C ) × 2(L y + 2L C ) × 2(L z + 2L C ) (with L i rounded up to the nearest multiple of L C ). At the end of the process, the desired sub-network in a domain of size L x × L y × L z was extracted from the center of the large domain.

S2 Results in lattice networks
Results for the two ordered, regular lattice networks introduced in Section 2.2.2 and shown in Figure 1Dare presented. These networks are referred to respectively as the CLN (with 6-connectivity) and the PLN (with 3-connectivity). After a study of the scaling properties of these networks, we present the results of metrics for chosen domain size and scaling.

S2.1 Scaling of lattice networks
We first present an analytical investigation into the scaling of these networks with domain size or with mean length L. Recall that in both lattice networks, the length of the elementary cubes, L C , is linearly proportional to L (L C = L in the CLN and L C = 3L in the PLN). The Voronoi-like synthetic networks, despite exhibiting greater randomness, also follow a similar semi-ordered structure controlled by the characteristic length L C . Thus we make the hypothesis that these metrics scale with L C in the same way, and use this to inform the scaling studies in Section 3.2.

S2.1.1 Morphometrical metrics
The mean vessel length was L for both lattice networks. The length density was 3L/L 3 for the CLN and 30L/(3L) 3 for the PLN. Obviously, none of these metrics depended on domain size, in contrast to the loop metrics studied below.

S2.1.2 Topological metrics
The scaling of topological metrics with domain size was derived analytically. Loop metrics computed for lattice networks generated for the range of domain sizes shown in Figure S2 agreed perfectly with the analytical expressions derived in this section.
For both networks, we suppose that i, j and k are the number of unit cells in the x, y, and z directions respectively. Both lattice networks had two modes of loops, α = 1, 2, with m α edges per loop. We then define the mean number of edges per loop, N edge/loop (i, j, k), as the total number of loop edges, divided by the total number of loops, i.e.: where N α loop (i, j, k) are the number of loops corresponding to mode α.
Similarly the mean loop length, L loop , was given by: and thus converged with domain size in the same way as N edge/loop (i, j, k), but also scaled with L. The mean number of loops per edge, N loop/edge (i, j, k), is the total number of loop edges divided by the number of interior (i.e. non-boundary) edges, N edge,int (i, j, k): These metrics were derived for the specific geometries of the CLNs and PLNs as follows.
Cubic lattice network The CLN had two modes of loops with 4 and 6 edges respectively (m = {4, 6}).
The number of loops with 4 edges, N 1 loop (i, j, k), followed: Similarly, the number of loops with 6 edges, N 2 loop (i, j, k), obeyed: where R(x) is the ramp function, defined as: Using Equation S1, the mean number of edges per loop for a cube of side length n = i = j = k and n ≥ 2 was given by: The mean number of edges per loop converged to within 5% of the converged value for networks of 4 3 unit cubes ( Figure S2a).
Convergence was slow: this metric was only converged to within 5% of this value for networks of 36 3 unit cubes or greater ( Figure S2b).
N edge/loop quickly converged to within 5% of this value for networks of 2 3 unit cubes ( Figure S2a).
N loop/edge converged more quickly for the PLN than the CLN, and was within 5% of its converged value for networks of 11 3 unit cubes or more ( Figure S2b).
Thus, in both lattice networks the number of loops per edge was especially sensitive to finite-size effects. Nonetheless, none of the loop metrics except the mean loop length depended on the choice of L (or equivalently L C ), and thus are purely topological metrics.

S2.1.3 Permeability
Cubic lattice network In the CLN, imposing a pressure gradient in the x-direction yields zero flow in vessels orientated parallel to the y-or z-axis: the network is reduced to an array of tubes of length L x

S4
where L x is the ROI size in the x-direction. The global flowrate Q x is: where µ is the effective viscosity, d is the diameter (which here is uniform), and ∆P is the pressure drop. N is the number of boundary vessels on the face x = 0, expressed in terms of L, the edge length: From Equation (1), this yields: Rearranging and using the definition for the domain area perpendicular to the x-axis, A x = L y L z , gives: Figure S3: Flow distribution in an elementary motif of the PLN with a pressure gradient from left to right. The magnitude of the flow in the green vessels was half of that in the red vessels, and zero in the dark blue vessels.
Periodic lattice network In the PLN, imposing a pressure gradient in the x-direction with uniform diameters also yields a simple flow distribution (Figure S3 and legend). The global flowrate obeys the same law in Equation S4, but N follows: In the minimal example shown in Figure S3, L y = 3L and L z = 3L giving N = 2. Substituting this expression into Equation (1) yields: Thus, to obtain the same permeability in both lattice networks, L in the PLN would have to be √ 2/3 = 0.47 times that of the CLN. Finally, for both lattice networks, the permeability scaled with d 4 and 1 /L 2 (or equivalently 1 /L 2 C )and did not depend on domain size (assuming uniform diameters).

S2.2 Lattice networks scaled to match mouse ROIs
In the CLN, L = 67µm was chosen to match the length density in mice. A network of 4 × 4 × 4 unit cubes was generated with dimensions 268 × 268 × 268µm 3 , to be as close as possible to the size of the mouse ROIs will containing whole unit cubes. In the PLN L = 41µm; with the unit cube side length being L C = 3L, a network of 2 × 2 × 2 unit cubes was generated with dimension 246 × 246 × 246µm 3 . Results in these networks are given in Table 1.
Space-filling metrics In the box-counting analysis, the cut-off lengths for both lattice networks were of a similar order of magnitude to the other networks, but in contrast, both had linear regimes at small scales. The mean EVD was very close to the mouse data for both lattice networks, while the maximum EVD was lowest in the CLN, and highest in the PLN. The mean local maxima of EVDs were 22% and 61% higher in the PLNs and CLNs respectively.

Morphometrical metrics
The PLN had a similar mean length to that in mice, while that in the CLN was almost 50% higher. By construction, the length distributions in these networks were drastically different to the mouse and synthetic networks, with only one mode for each network, and hence zero SD. Due to the choice of L as described above, the mean length density in both networks was very close to that in mice. As a result, in the PLN both the edge density and vertex density were close to those in mice, whereas the CLN had a much lower edge and vertex density. Both lattice networks had a much lower density of boundary vertices (36% and 62% respectively).
Topological metrics By construction the PLN had no multiply connected vertices, while all junctions in the CLN had 6 connections. The distri-S5 bution of loops was very different for the lattice networks, with fewer edges per loop (9 and 5.14 respectively) on average, with values restricted to two modes per network, as described above in Section S2.1.2. The mean loop length was also lower in both lattice networks than in mice.
Flow metrics The mean velocity was 45% higher in both the lattice networks. Similarly, the mean permeability in the PLNs and especially the CLNs was much higher than in mice (almost 30% and 120% higher respectively), and completely isotropic due to their regular construction.
Mass transfer metrics Unsurprisingly, the distribution of capillary transit times for lattice networks were very different to those in the synthetic and mouse ROIs ( Figure 8E). The exchange coefficient h was 26.5% and 32.9% higher in the PLNs and CLNs respectively compared mice.
Robustness metrics For the PLN, the histogram of pre-to post-occlusion flow ratios was very different to the distributions for the mice and synthetic networks, due to its highly organized architecture ( Figure 8F), and was restricted to 3 modes. The mean flow ratio was approximately 10% lower than in mice. The CLN was not included because none of its vertices were converging bifurcations.
In conclusion, while some simple morphometrical metrics in the lattice networks aligned well with those in mice, their highly organised structure meant that the distribution of metrics was completely different, and topological and functional metrics also differed greatly from the more randomized networks. Thus, such simple networks cannot be used to make accurate predictions of functional properties in the cerebral capillaries.  0.77 ± 0.03 0.76 ± 0.01 Post-occlusion flow ratio (diverging; branch A) 0.28 ± 0.05 0.28 ± 0.04 Table S1: The geometrical, topological and functional metrics calculated here, for human and synthetic with L C = 90µm ('S90') networks. Results are presented as mean ± S.D. over the N ROIs studied for each network type. Permeabilities, velocities and transit times were calculated with uniform diameters of 5µm.   Table S2: In the synthetic networks, results of a standard linear regression on a range of metrics as a function of L C to the appropriate power, i.e. following the expression: αL n C + β, where n is the power appropriate for each metric. The coefficient of determination, R 2 , was very good for all fits.