Network community structure and resilience to localized damage: application to brain microcirculation

In cerebrovascular networks, some vertices are more connected to each other than with the rest of the vasculature, defining a community structure. Here, we introduce a class of model networks built by rewiring Random Regular Graphs, which enables to reproduce this community structure and other topological properties of cerebrovascular networks. We use these model networks to study the global flow reduction induced by the removal of a single edge. We analytically show that this global flow reduction can be expressed as a function of the initial flow rate in the removed edge and of a topological quantity, both of which display probability distributions following Cauchy laws, i.e. with large tails. As a result, we show that the distribution of blood flow reductions is strongly influenced by the community structure. In particular, the probability of large flow reductions increases substantially when the community structure is stronger, weakening the network resilience to single capillary occlusions. We discuss the implications of these findings in the context of Alzheimer's Disease, in which the importance of vascular mechanisms, including capillary occlusions, is beginning to be uncovered.


I. INTRODUCTION
Cerebral hypoperfusion, i.e. the decrease of cerebral blood flow, is a common feature of many brain diseases, including neurodegenerative diseases, such as Alzheimer's Disease (AD) [1,2], and cerebrovascular diseases, such as hypoperfusion dementia [3]. Hypoperfusion is a key player in the onset and progression of cerebrovascular diseases [3,4] and has been considered until recently as a consequence of neurodegeneration in AD [5]. However, this view is now debated [1]. In human patients, cerebral blood flow indeed decreasesin a statistical and epidemiologic sense -before neurotoxic waste accumulate in the brain and before any measurable cognitive deficits [6]. Moreover, occlusions of capillary vessels by white blood cells (neutrophils) have been observed in animal models of AD before the accumulation of amyloid β , the main neurotoxic protein forming deposits (plaques) in AD brains [7]. Despite the small proportion of occluded vessels (from 1 to 4%), the pharmacological removal of these neutrophils led to a significant increase in blood flow and improved the cognitive performance of the animals. At a later stage , i.e. when the animals already showed plaques, extensive capillary constrictions have also been observed [8]. The exposure of pericytes, i.e. active mural cells wrapping around the capillaries, to increased concentrations of amyloid β has been shown to induce their contraction. This contributes to a positive feedback loop, where decreased cerebral blood flow not only triggers biological pathways leading to increased amyloid β production in the brain, but also directly impairs its elimination by the flowing blood. This results in increased amyloid β accumulation in the brain, increased pericyte contraction and further hypoperfusion [8]. In parallel, hypoperfusion also directly compromises the brain's energy supply, with deleterious neurological consequences.
A central question in this context is: to what extent could a small proportion of vessel occlusions trigger the above posi-tive feedback loop, contributing to AD onset and progression? Using highly resolved simulations of blood flow in anatomically realistic microvascular networks from human and mice, it was previously shown that, on average, cerebral blood flow decreases linearly with an increasing proportion of capillaries occluded at random, up to 20%, i.e. without any threshold effect [7]. Thus, on average, each single capillary occlusion has a similar, and cumulative, contribution to the blood flow decrease at the scale of the network. Here, we focus on the variability of this contribution and seek to determine how it is controlled by the fundamental topological properties of cerebrovascular networks. In particular, we investigate the role of network communities, i.e. substructures with vertices more connected to each other than to other vertices, which have been identified in such networks [9,10] and many other real-life networks [11,12]. In social networks, for example, these communities represent groups of users with specific affinities, the identification of which is an important area of research [13,14]. Such communities are also important to understand the organization of the World Wide Web [15,16], as they influence the various page ranking algorithms.
In order to study the role of communities in brain microvascular networks, we take a step back and adopt a theoretical point of view, using model networks of increasing complexity, as illustrated in Fig. 1, but with a single inlet and outlet, and vessels with identical unit conductivity. On the one hand, this abstract approach enables the network community structure to be controlled. On the other hand, it enables analytical expressions to be derived, or averages over many realizations of similar graphs with same properties to be performed, in order to analyze the variability of the blood flow reduction induced by the occlusion of a single vessel. We first consider simple ideal graphs known as Random Regular Graphs (RRGs) in network theory [17]. These graphs do not need to be embedded in the physical space and they are structureless. They nevertheless reproduce one of the main topological properties of cerebrovascular networks, in which most ver-tices have the same connectivity (or degree), equal to three (see e.g. [9,18,19]). Moreover, many realizations of RRGs of arbitrary size can be easily constructed, and they locally behave like trees, enabling analytical derivations that provide insight on their asymptotic behavior in the limit of large sizes. We then modify this ideal RRG model to provide a simple generation scheme that enables the strength of the communities to be controlled by rewiring together a finite number of elementary RRGs. As a third model, we use random networks constructed from Voronoi diagrams of sets of points homogeneously distributed in 3D space, following [19]. Such spatial networks are locally randomized but homogeneous at the network scale, and reproduce both the structure (morphology and topology) and function (flow, blood/tissue exchange and robustness to capillary occlusions) of brain capillary networks. Finally, we also consider the intracortical vascular network from the mouse parietal cortex (15,000 vessels in a 1mm 3 region) used in [7]. In these last two anatomically realistic networks, by contrast with [20,21], we neglect the contribution of vessel morphology, including distributions of diameters and lengths, and impose unit conductivity in all edges.
We use the above models to study the impact of single edge removal, i.e. equivalent to vessel occlusion. We show that the resulting flow reduction at network scale can be expressed as a function of the initial flow rate in the occluded vessel and a topological quantity, both of which display probability distributions with large tails. As a result, we show that the distribution of blood flow reductions may display unexpectedly large values, the probability of which increases when the community structure is stronger. Such results indicate that the topology of biological networks, including their community structure, is important to assess their functional properties, especially their biological resilience understood here as their ability to maintain functionality in the event of localized damage [22].
The paper is organized as follows. In Section II, we introduce the network models that will be considered, and investigate in Section III their topological properties compared to anatomically realistic networks. In Section IV, we consider blood flow through these networks. We uncover in particular the quantities controlling the distribution of flow reductions induced by the removal of a single edge and highlight how the network topological properties influence this distribution. Finally, in Section V, we discuss these findings and their implication for brain pathophysiology.

II. NETWORK MODELS
Network models have been used in many fields [23][24][25][26][27][28], where they proved useful to distinguish the main architectural properties of real-life complex networks that strongly impact their function from peculiarities which may depend on particular instances but do not change much the function.
Here, we start from the observation that, as stated in the Introduction, most vertices in cerebrovascular networks have three neighbours (either one vessel branches into two other vessels or two vessels merge into another one) [9,19]. The simplest graph model one can think of to study such networks is a graph model only reproducing this feature, i.e. with constant connectivity. Such graphs are known as regular graphs. Thus, the first model we consider is random regular graphs (RRGs), where each vertex has the same connectivity z, as illustrated in Fig. 1, left (see Appendix A for details).
RRG graphs by construction do not present specific substructures. However, in many real-life networks, including cerebrovascular networks, groups of vertices may have more links to each other than to other vertices [9,10], defining substructures, or communities. In order to account for such communities, we also consider a slightly more elaborate model, which we call rewired RRG model. It is obtained by generating n c independent RRG graphs, which may all have the same size or have a given heterogeneous size distribution, and rewiring pairs of edges at random, as illustrated in Fig. 1, middle (see Appendix A for details). In particular, we will denote by R k a subfamily of rewired RRGs built from a set of RRGs whose size is distributed to reproduce the communities of cerebrovascular networks, as further introduced in Section III B, and with k rewirings.
Both RRG and rewired RRG models are ideal graphs of infinite dimension [17,23,24,29], which are not embedded in the physical space (i.e. vertices have no a priori physical spatial coordinates). Thus, to account for the 3D structure of the brain, we consider a phenomenological model, constructed from edges of Voronoi cells in the three-dimensional space, which reproduce both the structure and function of brain capillary networks [19] (see Appendix A for details on the construction). The topology of a typical realization V of such graphs, that we call Voronoi graphs, is displayed in Fig. 1, right.
Finally, we also construct a graph based on experimental data from a mouse brain previously acquired by [9,18] (see Appendix A for details). By contrast to the previous graph, this experimental graph not only includes the space-filling capillary vessels, but also the tree-like penetrating arterioles and ascending venules [30]. Noteworthy, because these latter graphs are obtained from real 3D structures, vessels intersecting the edges of the considered region are cut. We recover the 3-connectivity of the graph by removing recursively all dangling vertices (vertices of connectivity one), and we keep only the largest connected component.

A. Distribution of loop lengths
The topology of a graph is often characterized by its mean connectivity [24]. However, graphs with mean identical connectivity may have very different topologies, from binary trees to periodic regular lattices or scale-free graphs. In the latter, hubs allow to connect any two vertices with a very short path [23], while in trees, there are no shortcuts. The distribution of loop length enables these different behaviors to be distinguished, while getting additional insight on network re-   dundancy. Moreover, it has been carefully studied in brain microvascular networks [9,19].
Here, following [19], we define a loop as the shortest path going from one vertex i to itself through two given neighbours of that vertex. For a vertex with z neighbours, there are z 2 loops. The loop length L is given by the total number of vertices (or edges) in the loop. In Fig. 2, we display the loop length distribution P(L) for the considered network models. As shown in the top panel, for a RRG of size N, the distribution is peaked around 2 ln N. This result is consistent with the RRGs being locally treelike [31,32], since in the limit N → ∞, the relative number of loops of fixed size p goes to 0 for any arbitrary integer p, making the graph effectively loopless. In the same way, the loop distribution of a collection of several independent elementary RRGs of size n 0 is peaked at 2 ln n 0 . Under increased rewiring, the distribution progressively shifts towards the one of a RRG of size N = n c n 0 . Thus, rewired RRGs continuously interpolate between a distribution of loops peaked at a value only depending on the substructure size n 0 to one peaked at a value only depending on the total network size N.
Interpreting the loop length distribution of the Voronoi and mouse networks is less trivial. As shown in the bottom panel of Fig. 2, orange and blue lines, respectively, they are peaked at a similar value, irrespective of their different network sizes, suggesting that, for these networks, the dominant factor is the community size and not the network size. Moreover, the mouse network has a slightly larger proportion of loops with lengths above the peak compared to the Voronoi graph V (orange line). This is consistent with earlier findings that the Voronoi graphs reproduce well the loop length distribution of cerebral capillary networks [19] and the additional contribution of tree-like arterioles and venules in the mouse network [9,33]. However, these distributions are larger than the typical distributions obtained in Fig. 2 top with rewired RRGs of equal size. Such distributions can be roughly reproduced by considering elementary RRGs with heterogeneous sizes (red lines in bottom panel). This suggests that the simple proposed model of RRG rewiring is sufficiently versatile to generate ideal graphs reproducing the topological properties of intracortical microvascular networks, by contrast to single RRGs, which do not account for the underlying substructures. In the next subsection, we thus focus on the community structure of these different networks.

B. Communities
As mentioned in the Introduction, communities are an ubiquitous feature of many real-life networks, and correspond to groups of vertices that are more likely to be connected together than with the rest of the graph. Such communities have been previously identified in brain microvascular networks [9,10]. To identify these communities, we maximize the modularity µ(C ) of all possible partitions C , following [11,12,34] and, in the case of intracortical networks, [9,10]. The quantity µ(C ) compares the probability of having an edge within a given subset of this partition with the probability expected by chance, i.e. from connections randomly chosen under the constraint of maintaining the graph connectivity z (see Appendix B for details). This quantity lies  between − 1 2 and 1, and is positive if the partition C has some relevance as a graph community structure. We call modularity µ of the graph the maximal modularity over all partitions, or equivalently the modularity of the optimal partition (see Appendix B). Noteworthily, even for graphs such as single RRGs without any built-in substructures, the optimization process finds the specific partition of the graph which maximizes the modularity. The value of µ(C ) for this particular partition is usually well above zero. Thus, for single RRGs, we obtain an average modularity µ RRG ∼ 0.661 when averaged over many realizations. This value corresponds to purely random fluctuations, which create random clusters of vertices more tightly bound together than with the rest of the graph. Higher values are expected for the graph community structure to have any relevance. In this case, we directly use the modularity to assess the strength of the obtained community structure. For the Voronoi graph we find µ V ∼ 0.827 and for the mouse graph µ mouse ∼ 0.846. Besides modularity, alternative quantities have also been used for that purpose, for example the exponent of the scaling law relating the number of inter-versus intra-community edges [9]. However, the range of exponents characterizing strong communities depends on the dimension of the space in which the graph is embedded [61], making it difficult to compare real-life networks and ideal graphs such as RRGs or rewired RRGs.
Below, we first examine the community structure of the Voronoi and mouse networks, and subsequently discuss how to incorporate these findings into the rewired RRG model. The community structures of the Voronoi and mouse networks are displayed in Fig. 3. Consistent with a high modularity value, the number of edges within communities is much larger than the number of edges connecting different communities. Moreover, these communities are highly segregated in space, i.e. with few overlap, as displayed in Fig. 4. In this figure, we use the spatial information about vertices of our 3D graphs to represent each community by a sphere, centered at the barycenter of its vertices and with a radius equal to the standard deviation of vertex positions around that barycenter. Not surprisingly, this suggests that in such systems, commu-  nities are a manifestation of the spatial organization of the network, and that vertices in a given spatial region make extensive connections with the neighboring vertices. By contrast, in pure RRGs, the community structure is mainly irrelevant. Communities in RRG graphs are highly connected to each other (Fig. 5, left) and the average modularity is much lower (µ RRG ∼ 0.661), even if, as said above, the maximization procedure always produces some spurious community structure by construction.
Relevant communities are recovered for the rewired RRGs (see e.g. Fig. 5, right, where communities associated to rewired RRGs of heterogeneous sizes are displayed). From the results presented in the previous Subsection, we expect that increasing the number of rewirings will result in approaching the behavior of a single, larger, RRG. This results in weaker communities, i.e. a network with lower modularity, as illustrated in Fig. 6, and an increased ratio of inter-versus intra-community edges. Below, we seek to estimate the number of rewirings needed to construct heterogeneous rewired RRGs with both community structures (defined by their number and size) and modularity matching those of a given reallife network.
For example, in the case of Voronoi V , the network has 1422 intra-community edges and 182 inter-community edges, and thus the ratio between inter and intra edges is ≈ 0.128. Choosing a rewired RRG with elementary RRG sizes given by the communities of V , we find a linear dependence of this ratio with the number of rewirings (inset of Fig. 6). To match the inter-over intra-edge ratio of V , this yields, via the fit of Fig. 6, n r = 74.79 ≈ 75 rewirings. The so obtained random graphs will be denoted by R 75 in the remaining of this paper. Their average modularity (obtained by averaging over 100 realizations) is equal to µ R 75 0.817 ± 0.0065, which is indeed very close to the modularity of the Voronoi and mouse graphs [62].
Noteworthy, the resulting distribution of loop lengths P(L) for such a number of rewirings matches well the distribution of loop lengths of Voronoi V (see Fig. 2), despite the fact that the 3D structure is entirely absent from the rewired RRG model. This similarity also manifests in the community structure displayed in Figs. 3 and 5.
The results of this Section show that topological features such as the community structure are important properties of intracortical microvascular networks. These properties can be implemented in the RRG model by adding a community structure and a certain amount of rewiring. Such models can describe correctly the loop distribution and the modularity of the real networks, while having no other features left in the model. They thus enable averages over all realizations of these networks, which all have the same structure, to be performed, a useful tool to get statistically significant quantities and assess the variability of the results. In fact, the rewired RRG ensemble corresponding to a given number of rewirings R k is the subset of the RRG ensemble which describes the universal behavior of regular graphs having a given heterogeneous com-munity size distribution and a given community strength. The Voronoi and mouse networks are specific instances of this ensemble and turn out to display topological quantities which are in line with the mean result of the equivalent rewired RRGs. We will use these models in the next Section to study the resilience of the above networks to single-edge occlusion and assess its dependence on modularity.

IV. BLOOD FLOW THROUGH THE NETWORK
We now connect the previous graphs with the outer world by adding a single inlet edge and a single outlet edge (see Fig. 7). Following [20][21][22], we also ignore the complex rheology of blood, so that imposing a constant pressure (potential) difference between the tips of these edges results in establishing a stationary flow field throughout the network, that can be obtained by inverting a linear system of equations [21,35]. In the present Section, we shall examine how occluding a single vessel affects the total flow rate transported through the network.

A. Definitions
Starting from a graph G with M = N − 2 vertices, we pick up two edges at random, and add two vertices I and O in the middle of these edges. To these vertices, which now are the N − 1 th and N th vertices of the graph, we respectively connect an inlet vertex I and an outlet vertex O . Note that, in that way, a graph with constant connectivity z = 3 keeps that property.
Let p i denote the potential (pressure) at vertex i for 1 i N. The local flux (flow rate) from i to j is defined as where γ is the matrix of conductivities.
where γ i j is the conductivity of edge i j. For i = N − 1 or N, there is an additional flux γ ii (p i − p i ), where p i is the corresponding imposed potential (p I or p 0 ), yielding where γ ii is the conductivity of the newly added edges, which will be denoted by γ I for the inlet and γ O for the outlet in the following. Let A be the N × N matrix defined by Laplacian matrix, and its spectral properties have been extensively studied [36,37]. We define vector b as the vector with entries −γ ii p i at position i = N − 1 and N and 0 elsewhere, that is, b = (0, . . . , 0, −γ I p I , −γ O p O ) T . Then, from Eqs. (2) and (3), the vector p is the unique solution of Ap = b. Generically, A is invertible (note that the sum over k in (4) includes i ), so that p is explicitly given by p = A −1 b.
In Fig. 7, we illustrate the fluxes q i j for one inlet point I on the left and one outlet point O on the right of a single RRG. By mass conservation, the total flow rate Q transported through the network is equal to the inlet and outlet fluxes. Thus, with the above notations, where Γ denotes the overall network conductance. We also denote by δ P = p I − p O the fixed value of the potential difference between the inlet and outlet. Then Q = Γδ P. We will now examine how Q, or equivalently Γ, is affected by the removal of one edge.

B. Flow reduction induced by removing one edge
We now denote by Q i j the total flow rate transported through the network when edge i − j is removed, and by δ Q i j = Q − Q i j the corresponding flow reduction. This flow reduction can be expressed as follows (see Appendix C) where q i j is the initial flux in the removed edge and In what follows, for simplicity, we will only consider the case where the pressure reduction δ P = p I − p O = 1 is fixed, so that Q = Γ, and where all conductivities γ i j are taken equal to 1. Equation (5) reduces to Equation (6) states that the total flow reduction induced by removing one edge i − j is quadratic in the initial flux q i j through this edge, and inversely proportional to t i j .  Noteworthy, the term t i j in Eq. (6) is analogous to expressions known as Line Outage Distribution Factors (LODF) describing the response to outage of an edge in an electric power grid [38]. Similar expressions have also been derived in the broader context of network theory [39], including for the study of single edge removal in space-filling model networks of biological vasculature [22].
Below, we further examine the properties of these two quantities, t i j and q i j , governing the flow reduction.

C. Cauchy laws for the quantities governing the flow reduction
It is first possible to get an insight into the denominator t i j in Eq. (6) by considering the asymptotic case of very large regular trees (see e.g. [40,41]). For that purpose, let us denote by G(x) = (A − x1) −1 the Green function of matrix A, so that A −1 = G(0). In the case where the off-diagonal elements of A are those of the adjacency matrix of an infinite regular tree, one can find a recursive relation between the diagonal elements of G(x). The latter expresses G ii (x) as a function of the G j j (x), with j the children of i, and of the diagonal elements of A (see e.g. [40,41]). For infinite trees, all vertices i are at the root of statistically identical trees. Therefore G ii (x) behaves as a random variable G and all G j j (x) behave as independent random variables G p with the same distribution as G. These are solutions of the recursive equation where z is the graph connectivity and c is a random variable distributed in the same way as the diagonal elements of A.
When the graph is no longer an infinite tree but a finite random regular graph, Eq. (7) is only approximate.
It is known in statistical physics [40,42] that the mean-field solution of this recursive equation amounts to approximating the exact solution (in the form of a probability distribution of G ii (x) as a random variable) by symmetric Cauchy distributions. As a result, (A −1 ) ii is expected to follow a Cauchy distribution, so that, by stability of Cauchy distributions, the quantity t i j = 1 + (A −1 ) ii + (A −1 ) j j − 2(A −1 ) i j , which corresponds to the denominator in (6), is also expected to follow a Cauchy distribution.
In Fig. 8, we plot the distributions of t i j for all our networks models, from RRGs to the mouse network. Whatever the considered graph, these distributions are indeed well-fitted by Cauchy distributions where t 0 is the median value and a is the half width at half maximum (HWHM) [43]. Besides, in all cases, we numerically find that the median value t 0 is close to 1/3, and thus the prefactor of q 2 i j in (6) is close to 3. This result can be recovered by generalizing the reasoning of [21,44], initially introduced to study the dilute regime of bond percolation. For regular networks of connectivity z with a single inlet and outlet, we show in Appendix D that the flow reduction can be approximated by in the limit of infinite size and under the assumption that the graph is isotropic at the inlet and outlet. The inset of Fig. 8 shows that this quadratic law is well verified whatever the considered graph. In our case, z = 3 and thus the prefactor of (p i − p j ) 2 in (9) is z z−2 = 3, recovering 1/3 as the median value t 0 of the Cauchy distribution for t i j .
In addition, Cauchy laws have fat tails. As a result, the occurrence probability of small and large values of quantities following Cauchy distributions is much larger than for normal distributions. For instance, the probability of values smaller than t 0 −3a, where t 0 denotes the median and a is the HWHM, is ≈ 0.0002 for a Gaussian distribution whereas it is ≈ 0.1 for a Cauchy distribution, i.e. five hundred times higher. Thus, there is a number of values of the t i j which are very small. When such an edge i − j is removed, this translates into a prefactor of the flow reduction q i j which is particularly high, inducing a significant drop in the total flow Q. As a result, the larger the Cauchy distribution in Fig. 8, the larger the fluctuations observed around Eq. (9) in the corresponding inset.
Importantly, the half width a of the Cauchy distribution strongly depends on the community structure of the underlying graph. In Fig. 9, we show that for graphs of comparable size, the width increases with the graph modularity, significantly increasing the probability of large flow reductions resulting from the removal of a single edge. Interestingly, while the Voronoi and mouse networks follow the same trend as the rewired RRGs R k , the width of their Cauchy distribution is larger than that obtained for rewired RRGs, up to an order of magnitude for the mouse network. It is unlikely that this increase can be due to the random dispersion of results in the rewired RRG ensemble. Rather, in the same way that the R k s, despite being specific realizations of graphs which belong to the RRG ensemble, display on average higher modularities than more probable realizations directly drawn from the RRG ensemble (compare modularity of red versus black symbols in Fig. 9), we may expect that other topological peculiarities of the community structure in the Voronoi and mouse networks (e.g. heterogeneous distribution of community size, increased probability of inter-links connecting neighbor communities in the 3D space) might also systematically bias the width of the Cauchy distribution. Additional biological data are probably needed to explore this point.
We now turn to the distribution of q i j , the numerator in (6). This distribution can be understood based on similar arguments as above. From (1), local fluxes are indeed obtained by linear combinations of potentials p = A −1 b, leading to a Cauchy distribution in the case of infinite trees. Besides, because of flow reversibility, its median value is expected to be zero. This is consistent with the blood flow distributions computed in realistic intracortical networks from mice [45]. These distributions have fat algebraic tails (power laws with exponent -2) characteristic of Cauchy distributions, which can alternately be interpreted, based on hydrodynamic arguments, as emerging from dipole-driven flows on random networks. Their properties have been investigated in [45], but the impact of the community structure has not been considered. Following the same approach as above, we find however that, for graphs of comparable size, modularity affects significantly less the width of the distribution of q i j compared to the one of t i j (see Fig. 10).

D. Distribution of flow reductions
After having determined the distribution of the denominators and numerators in Eq. (6), we now turn to the distribution of flow reductions δ Q i j , which controls the variability of the global network response to the removal of a single edge. There is no general result which describes this distribution for both q i j and t i j following Cauchy distributions. In Appendix E, however, assuming that the distribution of flow reductions is mainly controlled by the distribution of flow rates within the network, we show that for a regular graph and in the limit of large sizes, the left tail of the distribution of ln δ Q i j is expected to follow P(ln δ Q) ∼ 1/δ Q −1/2 and the right tail is expected to follow P(ln δ Q) ∼ 1/δ Q 1/2 .
As shown in Fig. 11 (insets), this asymptotic behavior is indeed correctly reproduced for the elementary RGGs (black lines), so that the distribution of ln δ Q is approximately symmetric. For these structureless graphs, it is also obvious that the position of its maximum, denoted by ln δ Q, strongly depends on the network size (compare black lines in top (N 1 =1176) versus bottom (N 2 =8720) panel in Fig. 11). As shown in Appendix E for large RRGs, we expect δ Q to scale as 1/N 2 . Thus, the two distributions obtained for RRGs with sizes N 1 = 1176 and N 2 = 8720 should match by a translation of the former by N 2 2 /N 2 1 to the left. This is verified in the bottom panel of Fig. 11 (superimposed black line and black dots) and suggests that, for these structureless graphs, the distribution of flow reductions is indeed controlled by the distribution of flow rates throughout the network. However, whatever the graph size, the distributions of ln δ Q become distorted towards larger values when the modularity increases, i.e. for graphs with substructures. As a result, the corresponding dis-   Figure 12: Probability of flow reductions larger than flow threshold δ Q 0 as a function of modularity. δ Q 0 is chosen to correspond to the third quartile of the distribution obtained for the RRG of same size, as displayed in Fig. 11 by dashed vertical lines. Same color code as Fig. 11 (RRG with N 1 = 1176 is the black circle, RRG with N 2 = 8720 is the black triangle. tributions slowly depart from the above predicted scalings, especially for large values of δ Q i j , and the location of their maxima increases. This is consistent with the concomitant in-crease of a, the width of the distribution of the denominator in Eq. (6), demonstrated in the previous Section.
To further quantify the effect of modularity, we plot in Fig. 12 the probability of flow reductions larger than a given threshold δ Q 0 , chosen to correspond to the third quartile of the distribution obtained for the RRG of same size [63]. From this arbitrary choice, this proportion is obviously 25% for both RRGs. It increases almost linearly for increasing modularities, up to ∼30% for the rewired RRGs. It even increases by almost twofold for the Voronoi and mouse networks. Consistent with the results presented in Fig. 9, this larger increase is likely due to topological peculiarities of the community structure in the Voronoi and mouse networks (e.g. heterogeneous distribution of community size, increased probability of interlinks connecting neighbor communities in the 3D space), but additional biological data are needed to explore this point.
Overall, in general, a stronger community structure is clearly associated with a higher probability of large flow reductions when removing one edge. As a result, in brain microvascular networks, the removal (or equivalently the occlusion) of one single vessel yields probabilities of large blood flow reductions that may be more than twice the prediction for graphs without substructures. Our interpretation of this result is that the community structure implies that there are subparts of the network less coupled with each other. If the occlusion affects a vessel linking two different communities, the flow will be very much affected since such inter-community edges are in a smaller number and the alternative flow paths can be very different.

V. CONCLUSION AND PERSPECTIVES
We have investigated the topological structure of cerebrovascular networks, based on an anatomical network extracted from a mouse brain, model capillary networks derived from space-filling Voronoi tessellations and different graph models belonging to the class of Random Regular Graphs. We have shown that the anatomical network contains substructures corresponding to communities of vessels, which are geographically localized in 3D. The proposed models, including such communities, have been obtained by rewiring RRGs. They thus belong to graph ensembles describing the universal behavior of regular graphs having a given community strength, and it is possible to tune the number of rewirings so as to obtain an average modularity in line with the modularity of any vascular network. Moreover, they provide a reasonable representation of the topology of cerebrovascular networks, as highlighted by their similar loop distributions. Noteworthily, whether such a generation scheme could be generalized to describe the topological properties of other vascular networks, especially pathological (e.g. tumor) ones, remains an interesting open question.
Using these simple network models, we have then studied how the strength of their communities affects the distribution of global flow reductions when an edge is removed, or, equivalently, when a single vessel is occluded. For that purpose, we have shown that the resulting flow reduction at net-work scale can be expressed as a function of the initial flow rate in the occluded vessel and a topological quantity. Both of them display probability distributions with large tails, resulting in a large probability of rare events. The presence of community structures even enlarges the tail of the distribution of the topological quantity, resulting in general in a distribution of flow reductions shifted towards larger and larger values when the community structure becomes more pronounced. In other words, the community structure weakens the network resilience to capillary occlusions, by increasing the probability of larger flow reductions, inducing a large variability of the impact of a single vessel occlusion depending on its location in the network.
The proposed theoretical approach neglects the heterogeneity of vessel conductivities. In other words, by contrast to [20][21][22]46], we ignored the contribution of vessel morphology, including distributions of diameters and lengths. We also ignored the complex rheology of blood [35,47,48]. Our results are nevertheless consistent with recent work modeling the impact of single capillary occlusion in highly resolved numerical simulations accounting for blood microvascular rheology in anatomically realistic microvascular networks [46], focusing on local flow reorganizations in the vicinity of the occluded vessel. Not surprisingly, the results, obtained for a total of 96 occluded capillaries representing less than 1% of capillaries, exhibited considerable numerical dispersion, which has been pointed out by the authors as a methodological difficulty. Nevertheless, the median volume of the region with flow reductions above 20% has been shown to increase by a factor 2.5 between capillaries with low initial flow rate and those with high initial flow rate. This suggests that the present theoretical framework, which enables the use of graph models focusing on network topology, keeps enough physical ingredients to be relevant.
While this does not sufficiently reduce the complexity of the problem to yield a complete theory relating the community structure to the shape of the distribution of flow reductions (providing e.g full analytical derivations), it still enables relevant asymptotic scalings to be deduced for all quantities controlling this distribution, which considerably helps to interpret the numerical results obtained in the mouse and Voronoi networks. Noteworthy, the last decades have seen tremendous progress of in vivo experimental techniques, including multiphoton microscopy and laser-based techniques that enable microvessels in rodent brains to be selectively occluded [49]. This offers large possibilities of data collection, which may be useful to validate our findings. In a complementary way, the present theoretical framework may help enrich data interpretation, e.g. by considering the impact of network communities on spatial flow redistribution [22] or by enabling the expected broad distribution of flow reductions to be taken into account in the statistical design of the studies.
Moreover, the present work provides a theoretical basis for future studies about the impact of multiple capillary occlusions [7] on cerebrovascular function. As mentioned in the Introduction, this may help to understand the interplay between hypoperfusion and amyloid-induced neurodegeneration in the onset and progression of AD [1]. We may in particular speculate that the microvascular community structure evolves with disease progression. For example, capillary occlusions at early stages of the disease, i.e. in healthy networks, may strengthen their community structure. This would increase the probability of larger blood flow reductions induced by further occlusions, providing an additional self-amplificatory mechanism in the positive feedback loop linking hypoperfusion and neurodegeneration in AD [1]. In the same way, different network organizations in different brain areas (e.g. primary versus secondary cortex, subcortical regions, hippocampus), which are being uncovered thanks to whole brain post-mortem vascular network reconstructions in rodents [50], may be a clue to explain their different vulnerabilities, e.g. understand why the hippocampus is one of the first damaged brain regions exhibiting cognitive deficits in AD. Long-term vascular remodelling, including capillary rarefaction, in normal or pathological aging [51] may also contribute to modify the community structure of vascular networks in the brain, thus providing an additional mechanism which may explain the considerable overlap between vascular pathology and AD [3,5].
To investigate the above assumptions, new datasets finely mapping the whole-brain vascular architecture of rodents in normal aging and at different stages of various brain diseases, including AD, are needed. Moreover, the present framework should be enriched to account for multiple occlusions, e.g. by introducing a perturbation approach valid in the dilute limit, where the removal of edges can still be considered independent of each other. Multiple network inlets and outlets should also be considered. The vascular network within the brain cortex is indeed fed and drained by a large number of penetrating arterioles and ascending venules [9,52], which may contribute to enhance the network resilience to capillary occlusions [22]. By contrast, occlusions of penetrating arterioles induce dramatic damage, as shown by a comprehensive series of in vivo experiments [9,30,[53][54][55]. This has been interpreted as resulting from insufficient compensatory collateral flow from other network inlets through the capillary bed, as reviewed in [52]. Our result suggest that it may alternately be understood as resulting from the hierarchical organization of the network, inlet (and outlet) vessels being those carrying the largest flow rates, thus leading to the largest flow reductions, both at network scale and, by extension, recursively in the neighborhood of the occluded vessel. This phenomenon is likely to be increased if conductance heterogenities are taken into account, such vessels displaying the largest conductances, leading to correlations between high flow and high conductance [57], while such correlations are negligible in the capillary network. In particular, studying whether these high conductance / high flow vessels have any relationship with the community structure would be an interesting extension of this work. Besides AD, this would open perspectives in the context of ischemic stroke, where neutrophil occlusions of up to 30% of capillary vessel have been recently discovered, preventing reperfusion after recanalization of the upstream cerebral artery [56].
Voronoi graphs. In addition to having mostly 3-connected vertices, the network of capillary vessels, i.e. the smallest vessels within the brain cortex, is space-filling [19,33]. This last property can be reproduced by constructing 3D Voronoi diagrams from sets of seed points randomly distributed under the strong constraint that there is only one seed point in each cube of a 3D cubic grid. However, these Voronoi diagrams have high connectivity, with many vertices of degree up to 5. By randomly merging, pruning or adding vertices following the geometrical constraints described in [19], we get 3D model networks statistically reproducing most of the morphological, topological and transport properties of brain capillary networks [19]. Because these networks are generated in a 3D cubic region, all boundary edges are cut and dangling, so that their outer boundary vertex is only 1-connected. In the present paper, we recursively remove these dangling edges, so that all remaining vertices are at least of degree 2. For simplicity, the resulting graphs are described as Voronoi graphs in the present paper. One of them, which we denote V , serves as an illustration throughout the paper.
Mouse graph. We use the graph description of a large postmortem dataset (∼ 1 mm 3 and ∼ 15, 000 vessel segments) from the mouse vibrissa primary sensory (vS1) cortex previously obtained by [9,18]. While this dataset contains vessel diameters and labels classifying vessels in arterioles, capillaries and venules, we discard this information and consider that all edges are equivalent, with unit conductivities. As above, we recursively remove all dangling edges.

Appendix B: Communities and modularity
A graph G with N vertices can be partitioned into communities, which are subsets C k of vertices. Intuitively, a community C k is such that vertices in C k are highly connected to one another, with comparatively fewer edges connecting them to vertices outside the community. There is however no unique answer to the question of what a meaningful partition of a graph is. Under reasonable assumptions, it is possible to construct many functions that quantify the relevance of a given partition for community detection [12].
One such function which has been very much used for that purpose in network theory is the modularity [34], which for a given partition C into subsets C k is defined as In this expression a i j is equal to 1 if there is an edge connecting vertices i and j, and to 0 otherwise, d i is the number of outgoing edges of i, m is the total number of edges, and the sum runs over all subsets C k in the partition. The sum ∑ i, j∈C k a i j gives (twice) the number of edges within the set C k . The number d i d j is the number of edges that could connect i and j if they were taken at random. The modularity thus compares the mean probability of having an edge of the graph G within a set C k to the probability of having such an edge in a graph G where all vertices have the same degree as in G but edges are chosen at random. The various possible functions characterizing the relevance of the partition lead to different methods for community detection in graphs (for a review see [12]). We used modularitybased clustering algorithms, which go through the space of possible partitions trying to maximize the modularity. In practice, we use the FindGraphCommunities function of Mathematica. The algorithm finds the optimal partition with the largest modularity for a given graph. This enables definition of the modularity of the graph, which we identify with the modularity of this optimal partition, or equivalently the maximal modularity over all partitions C : We then define the community structure of a graph as the partition with this maximal modularity (B2).
Appendix C: Removing one edge We recall that potentials are solution of Ap = b, where A is the matrix of conductances and b is the given by the righthand side of Eq. (2), namely b = (0, . . . , 0, −γ I p I , −γ O p O ) T (i.e. potentials p I and p O are imposed at inlet vertex I and outlet vertex O , respectively, see Fig. 7). We now examine the consequence, on the total flow rate Q, of removing one edge from the network. Let us first suppose, without loss of generality, that there is some edge connecting vertices 1 and 2 and that we remove it. Let us denote by A the matrix Eq. (4) with edge 1-2 removed. The new potentials p j in the absence of edge 1-2 are solution of A p = b.
Matrices A and A only differ by their upper-left 2 × 2 corner: the off-diagonal elements are indeed given by A 12 = γ 12 when edge 1-2 is present, and A 12 = 0 when it has been removed, while the diagonal elements change from A ii to A ii = A ii + γ 12 . Introducing the column vector u defined by u T = {1, −1, 0, ..., 0}, this can be reexpressed as A = A + γ 12 uu T .
Multiplying both members of this equation by vector b T on the left and b on the right, we get b T p = b T p − γ 12 p T uu T p 1 + γ 12 u T A −1 u , where p = (A + γ 12 uu T ) −1 b is the solution to the flow equation with edge 1-2 removed. From the definition of u, the scalar product of p and u is p T u = p 1 − p 2 . As for the scalar product b T p, we use the fact that Q = γ I (p I − p I ) = γ O (p O − p O ) = Γδ P, which leads to This equation is exact, and since there is nothing special about vertices 1 and 2, it remains valid for any arbitrary edge removed. Thus, in general, we have where δ Q i j = Q − Q i j , and t i j = 1 + γ i j [(A −1 ) ii + (A −1 ) j j − 2(A −1 ) i j ] only depends on the network topology and is dimensionless. This equation is homogeneous and leads to Eq. (5).
with a some constant. Besides, the root of the derivative of Eq. (E2) yields the maximum of the distribution p for δ Q = 3Q 2 c . Finally, in the limit of large sizes, RRGs of connectivity 3 with one additional inlet and outlet behave like the union of two balanced binary trees of equal height H [43], with roots corresponding to the inlet and outlet. Thus, we have N = 2(2 H − 1) − N l , N l = 2 (H−1) corresponding to the number of leaves that merge to connect the two trees. For such a graph, the distribution of absolute flow rate is fully described by the power-law regime, so that Q c is equal to the lowest value of the distribution, i.e. in the graph leaves. As a result, Q c = Q/N l . Combining these two equations leads to Q c = 3Q/(N + 1), so that, in the limit of large sizes, δ Q scales as 1/N 2 .